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   nlhdlr_bilinear.c
26   	 * @ingroup DEFPLUGINS_NLHDLR
27   	 * @brief  bilinear nonlinear handler
28   	 * @author Benjamin Mueller
29   	 */
30   	
31   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
32   	
33   	#include <string.h>
34   	
35   	#include "scip/nlhdlr_bilinear.h"
36   	#include "scip/cons_nonlinear.h"
37   	#include "scip/expr_product.h"
38   	#include "scip/expr_var.h"
39   	
40   	/* fundamental nonlinear handler properties */
41   	#define NLHDLR_NAME               "bilinear"
42   	#define NLHDLR_DESC               "bilinear handler for expressions"
43   	#define NLHDLR_DETECTPRIORITY     -10  /**< it is important that the nlhdlr runs after the default nldhlr */
44   	#define NLHDLR_ENFOPRIORITY       -10
45   	
46   	#define MIN_INTERIORITY           0.01 /**< minimum interiority for a reference point for applying separation */
47   	#define MIN_ABSBOUNDSIZE          0.1  /**< minimum size of variable bounds for applying separation */
48   	
49   	/* properties of the bilinear nlhdlr statistics table */
50   	#define TABLE_NAME_BILINEAR                 "nlhdlr_bilinear"
51   	#define TABLE_DESC_BILINEAR                 "bilinear nlhdlr statistics table"
52   	#define TABLE_POSITION_BILINEAR             14800                  /**< the position of the statistics table */
53   	#define TABLE_EARLIEST_STAGE_BILINEAR       SCIP_STAGE_INITSOLVE   /**< output of the statistics table is only printed from this stage onwards */
54   	
55   	
56   	/*
57   	 * Data structures
58   	 */
59   	
60   	/** nonlinear handler expression data */
61   	struct SCIP_NlhdlrExprData
62   	{
63   	   SCIP_Real             underineqs[6];      /**< inequalities for underestimation */
64   	   int                   nunderineqs;        /**< total number of inequalities for underestimation */
65   	   SCIP_Real             overineqs[6];       /**< inequalities for overestimation */
66   	   int                   noverineqs;         /**< total number of inequalities for overestimation */
67   	   SCIP_Longint          lastnodeid;         /**< id of the last node that has been used for separation */
68   	   int                   nseparoundslastnode; /**< number of separation calls of the last node */
69   	};
70   	
71   	/** nonlinear handler data */
72   	struct SCIP_NlhdlrData
73   	{
74   	   SCIP_EXPR**           exprs;             /**< expressions that have been detected by the nlhdlr */
75   	   int                   nexprs;            /**< total number of expression that have been detected */
76   	   int                   exprsize;          /**< size of exprs array */
77   	   SCIP_HASHMAP*         exprmap;           /**< hashmap to store the position of each expression in the exprs array */
78   	
79   	   /* parameter */
80   	   SCIP_Bool             useinteval;        /**< whether to use the interval evaluation callback of the nlhdlr */
81   	   SCIP_Bool             usereverseprop;    /**< whether to use the reverse propagation callback of the nlhdlr */
82   	   int                   maxseparoundsroot; /**< maximum number of separation rounds in the root node */
83   	   int                   maxseparounds;     /**< maximum number of separation rounds in a local node */
84   	   int                   maxsepadepth;      /**< maximum depth to apply separation */
85   	};
86   	
87   	/*
88   	 * Local methods
89   	 */
90   	
91   	/** helper function to compute the violation of an inequality of the form xcoef * x <= ycoef * y + constant for two
92   	 *  corner points of the domain [lbx,ubx] x [lby,uby]
93   	 */
94   	static
95   	void getIneqViol(
96   	   SCIP_VAR*             x,                  /**< first variable */
97   	   SCIP_VAR*             y,                  /**< second variable */
98   	   SCIP_Real             xcoef,              /**< x-coefficient */
99   	   SCIP_Real             ycoef,              /**< y-coefficient */
100  	   SCIP_Real             constant,           /**< constant */
101  	   SCIP_Real*            viol1,              /**< buffer to store the violation of the first corner point */
102  	   SCIP_Real*            viol2               /**< buffer to store the violation of the second corner point */
103  	   )
104  	{
105  	   SCIP_Real norm;
106  	   assert(viol1 != NULL);
107  	   assert(viol2 != NULL);
108  	
109  	   norm = SQRT(SQR(xcoef) + SQR(ycoef));
110  	
111  	   /* inequality can be used for underestimating xy if and only if xcoef * ycoef > 0 */
112  	   if( xcoef * ycoef >= 0 )
113  	   {
114  	      /* violation for top-left and bottom-right corner */
115  	      *viol1 = MAX(0, (xcoef * SCIPvarGetLbLocal(x)  - ycoef * SCIPvarGetUbLocal(y) - constant) / norm); /*lint !e666*/ /*lint !e661*/
116  	      *viol2 = MAX(0, (xcoef * SCIPvarGetUbLocal(x)  - ycoef * SCIPvarGetLbLocal(y) - constant) / norm); /*lint !e666*/ /*lint !e661*/
117  	   }
118  	   else
119  	   {
120  	      /* violation for top-right and bottom-left corner */
121  	      *viol1 = MAX(0, (xcoef * SCIPvarGetUbLocal(x)  - ycoef * SCIPvarGetUbLocal(y) - constant) / norm); /*lint !e666*/ /*lint !e661*/
122  	      *viol2 = MAX(0, (xcoef * SCIPvarGetLbLocal(x)  - ycoef * SCIPvarGetLbLocal(y) - constant) / norm); /*lint !e666*/ /*lint !e661*/
123  	   }
124  	}
125  	
126  	/** auxiliary function to decide whether to use inequalities for a strong relaxation of bilinear terms or not */
127  	static
128  	SCIP_Bool useBilinIneqs(
129  	   SCIP*                 scip,               /**< SCIP data structure */
130  	   SCIP_VAR*             x,                  /**< x variable */
131  	   SCIP_VAR*             y,                  /**< y variable */
132  	   SCIP_Real             refx,               /**< reference point for x */
133  	   SCIP_Real             refy                /**< reference point for y */
134  	   )
135  	{
136  	   SCIP_Real lbx;
137  	   SCIP_Real ubx;
138  	   SCIP_Real lby;
139  	   SCIP_Real uby;
140  	   SCIP_Real interiorityx;
141  	   SCIP_Real interiorityy;
142  	   SCIP_Real interiority;
143  	
144  	   assert(x != NULL);
145  	   assert(y != NULL);
146  	   assert(x != y);
147  	
148  	   /* get variable bounds */
149  	   lbx = SCIPvarGetLbLocal(x);
150  	   ubx = SCIPvarGetUbLocal(x);
151  	   lby = SCIPvarGetLbLocal(y);
152  	   uby = SCIPvarGetUbLocal(y);
153  	
154  	   /* compute interiority */
155  	   interiorityx = MIN(refx-lbx, ubx-refx) / MAX(ubx-lbx, SCIPepsilon(scip)); /*lint !e666*/
156  	   interiorityy = MIN(refy-lby, uby-refy) / MAX(uby-lby, SCIPepsilon(scip)); /*lint !e666*/
157  	   interiority = 2.0*MIN(interiorityx, interiorityy);
158  	
159  	   return ubx - lbx >= MIN_ABSBOUNDSIZE && uby - lby >= MIN_ABSBOUNDSIZE && interiority >= MIN_INTERIORITY;
160  	}
161  	
162  	/** helper function to update the best relaxation for a bilinear term when using valid linear inequalities */
163  	static
164  	void updateBilinearRelaxation(
165  	   SCIP*                 scip,               /**< SCIP data structure */
166  	   SCIP_VAR* RESTRICT    x,                  /**< first variable */
167  	   SCIP_VAR* RESTRICT    y,                  /**< second variable */
168  	   SCIP_Real             bilincoef,          /**< coefficient of the bilinear term */
169  	   SCIP_SIDETYPE         violside,           /**< side of quadratic constraint that is violated */
170  	   SCIP_Real             refx,               /**< reference point for the x variable */
171  	   SCIP_Real             refy,               /**< reference point for the y variable */
172  	   SCIP_Real* RESTRICT   ineqs,              /**< coefficients of each linear inequality; stored as triple (xcoef,ycoef,constant) */
173  	   int                   nineqs,             /**< total number of inequalities */
174  	   SCIP_Real             mccormickval,       /**< value of the McCormick relaxation at the reference point */
175  	   SCIP_Real* RESTRICT   bestcoefx,          /**< pointer to update the x coefficient */
176  	   SCIP_Real* RESTRICT   bestcoefy,          /**< pointer to update the y coefficient */
177  	   SCIP_Real* RESTRICT   bestconst,          /**< pointer to update the constant */
178  	   SCIP_Real* RESTRICT   bestval,            /**< value of the best relaxation that have been found so far */
179  	   SCIP_Bool*            success             /**< buffer to store whether we found a better relaxation */
180  	   )
181  	{
182  	   SCIP_Real constshift[2] = {0.0, 0.0};
183  	   SCIP_Real constant;
184  	   SCIP_Real xcoef;
185  	   SCIP_Real ycoef;
186  	   SCIP_Real lbx;
187  	   SCIP_Real ubx;
188  	   SCIP_Real lby;
189  	   SCIP_Real uby;
190  	   SCIP_Bool update;
191  	   SCIP_Bool overestimate;
192  	   int i;
193  	
194  	   assert(x != y);
195  	   assert(!SCIPisZero(scip, bilincoef));
196  	   assert(nineqs >= 0 && nineqs <= 2);
197  	   assert(bestcoefx != NULL);
198  	   assert(bestcoefy != NULL);
199  	   assert(bestconst != NULL);
200  	   assert(bestval != NULL);
201  	
202  	   /* no inequalities available */
203  	   if( nineqs == 0 )
204  	      return;
205  	   assert(ineqs != NULL);
206  	
207  	   lbx = SCIPvarGetLbLocal(x);
208  	   ubx = SCIPvarGetUbLocal(x);
209  	   lby = SCIPvarGetLbLocal(y);
210  	   uby = SCIPvarGetUbLocal(y);
211  	   overestimate = (violside == SCIP_SIDETYPE_LEFT);
212  	
213  	   /* check cases for which we can't compute a tighter relaxation */
214  	   if( SCIPisFeasLE(scip, refx, lbx) || SCIPisFeasGE(scip, refx, ubx)
215  	      || SCIPisFeasLE(scip, refy, lby) || SCIPisFeasGE(scip, refy, uby) )
216  	      return;
217  	
218  	   /* due to the feasibility tolerances of the LP and NLP solver, it might possible that the reference point is
219  	    * violating the linear inequalities; to ensure that we compute a valid underestimate, we relax the linear
220  	    * inequality by changing its constant part
221  	    */
222  	   for( i = 0; i < nineqs; ++i )
223  	   {
224  	      constshift[i] = MAX(0.0, ineqs[3*i] * refx - ineqs[3*i+1] * refy - ineqs[3*i+2]);
225  	      SCIPdebugMsg(scip, "constant shift of inequality %d = %.16f\n", i, constshift[i]);
226  	   }
227  	
228  	   /* try to use both inequalities */
229  	   if( nineqs == 2 )
230  	   {
231  	      SCIPcomputeBilinEnvelope2(scip, bilincoef, lbx, ubx, refx, lby, uby, refy, overestimate, ineqs[0], ineqs[1],
232  	         ineqs[2] + constshift[0], ineqs[3], ineqs[4], ineqs[5] + constshift[1], &xcoef, &ycoef, &constant, &update);
233  	
234  	      if( update )
235  	      {
236  	         SCIP_Real val = xcoef * refx + ycoef * refy + constant;
237  	         SCIP_Real relimpr = 1.0 - (REALABS(val - bilincoef * refx * refy) + 1e-4) / (REALABS(*bestval - bilincoef * refx * refy) + 1e-4);
238  	         SCIP_Real absimpr = REALABS(val - (*bestval));
239  	
240  	         /* update relaxation if possible */
241  	         if( relimpr > 0.05 && absimpr > 1e-3 && ((overestimate && SCIPisRelLT(scip, val, *bestval))
242  	             || (!overestimate && SCIPisRelGT(scip, val, *bestval))) )
243  	         {
244  	            *bestcoefx = xcoef;
245  	            *bestcoefy = ycoef;
246  	            *bestconst = constant;
247  	            *bestval = val;
248  	            *success = TRUE;
249  	         }
250  	      }
251  	   }
252  	
253  	   /* use inequalities individually */
254  	   for( i = 0; i < nineqs; ++i )
255  	   {
256  	      SCIPcomputeBilinEnvelope1(scip, bilincoef, lbx, ubx, refx, lby, uby, refy, overestimate, ineqs[3*i], ineqs[3*i+1],
257  	         ineqs[3*i+2] + constshift[i], &xcoef, &ycoef, &constant, &update);
258  	
259  	      if( update )
260  	      {
261  	         SCIP_Real val = xcoef * refx + ycoef * refy + constant;
262  	         SCIP_Real relimpr = 1.0 - (REALABS(val - bilincoef * refx * refy) + 1e-4)
263  	               / (REALABS(mccormickval - bilincoef * refx * refy) + 1e-4);
264  	         SCIP_Real absimpr = REALABS(val - (*bestval));
265  	
266  	         /* update relaxation if possible */
267  	         if( relimpr > 0.05 && absimpr > 1e-3 && ((overestimate && SCIPisRelLT(scip, val, *bestval))
268  	             || (!overestimate && SCIPisRelGT(scip, val, *bestval))) )
269  	         {
270  	            *bestcoefx = xcoef;
271  	            *bestcoefy = ycoef;
272  	            *bestconst = constant;
273  	            *bestval = val;
274  	            *success = TRUE;
275  	         }
276  	      }
277  	   }
278  	}
279  	
280  	/** helper function to determine whether a given point satisfy given inequalities */
281  	static
282  	SCIP_Bool isPointFeasible(
283  	   SCIP*                 scip,               /**< SCIP data structure */
284  	   SCIP_Real             x,                  /**< x-coordinate */
285  	   SCIP_Real             y,                  /**< y-coordinate */
286  	   SCIP_Real             lbx,                /**< lower bound of x */
287  	   SCIP_Real             ubx,                /**< upper bound of x */
288  	   SCIP_Real             lby,                /**< lower bound of y */
289  	   SCIP_Real             uby,                /**< upper bound of y */
290  	   SCIP_Real*            ineqs,              /**< inequalities of the form coefx x <= coefy y + constant */
291  	   int                   nineqs              /**< total number of inequalities */
292  	   )
293  	{
294  	   int i;
295  	
296  	   assert(ineqs != NULL);
297  	   assert(nineqs > 0);
298  	
299  	   /* check whether point satisfies the bounds */
300  	   if( SCIPisLT(scip, x, lbx) || SCIPisGT(scip, x, ubx)
301  	      || SCIPisLT(scip, y, lby) || SCIPisGT(scip, y, uby) )
302  	      return FALSE;
303  	
304  	   /* check whether point satisfy the linear inequalities */
305  	   for( i = 0; i < nineqs; ++i )
306  	   {
307  	      SCIP_Real coefx = ineqs[3*i];
308  	      SCIP_Real coefy = ineqs[3*i+1];
309  	      SCIP_Real constant = ineqs[3*i+2];
310  	
311  	      /* TODO check with an absolute comparison? */
312  	      if( SCIPisGT(scip, coefx*x - coefy*y - constant, 0.0) )
313  	         return FALSE;
314  	   }
315  	
316  	   return TRUE;
317  	}
318  	
319  	/** helper function for computing all vertices of the polytope described by the linear inequalities and the local
320  	 *  extrema of the bilinear term along each inequality
321  	 *
322  	 * @note there are at most 22 points where the min/max can be achieved (given that there are at most 4 inequalities)
323  	 *          - corners of [lbx,ubx]x[lby,uby] (4)
324  	 *          - two intersection points for each inequality with the box (8)
325  	 *          - global maximum / minimum on each inequality (4)
326  	 *          - intersection between two inequalities (6)
327  	 */
328  	static
329  	void getFeasiblePointsBilinear(
330  	   SCIP*                 scip,               /**< SCIP data structure */
331  	   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler, if levelset == TRUE, otherwise can be NULL */
332  	   SCIP_EXPR*            expr,               /**< product expression */
333  	   SCIP_INTERVAL         exprbounds,         /**< bounds on product expression, only used if levelset == TRUE */
334  	   SCIP_Real*            underineqs,         /**< inequalities for underestimation */
335  	   int                   nunderineqs,        /**< total number of inequalities for underestimation */
336  	   SCIP_Real*            overineqs,          /**< inequalities for overestimation */
337  	   int                   noverineqs,         /**< total number of inequalities for overestimation */
338  	   SCIP_Bool             levelset,           /**< should the level set be considered? */
339  	   SCIP_Real*            xs,                 /**< array to store x-coordinates of computed points */
340  	   SCIP_Real*            ys,                 /**< array to store y-coordinates of computed points */
341  	   int*                  npoints             /**< buffer to store the total number of computed points */
342  	   )
343  	{
344  	   SCIP_EXPR* child1;
345  	   SCIP_EXPR* child2;
346  	   SCIP_Real ineqs[12];
347  	   SCIP_INTERVAL boundsx;
348  	   SCIP_INTERVAL boundsy;
349  	   SCIP_Real lbx;
350  	   SCIP_Real ubx;
351  	   SCIP_Real lby;
352  	   SCIP_Real uby;
353  	   int nineqs = 0;
354  	   int i;
355  	
356  	   assert(scip != NULL);
357  	   assert(conshdlr != NULL || !levelset);
358  	   assert(expr != NULL);
359  	   assert(xs != NULL);
360  	   assert(ys != NULL);
361  	   assert(SCIPexprGetNChildren(expr) == 2);
362  	   assert(noverineqs + nunderineqs > 0);
363  	   assert(noverineqs + nunderineqs <= 4);
364  	
365  	   *npoints = 0;
366  	
367  	   /* collect inequalities */
368  	   for( i = 0; i < noverineqs; ++i )
369  	   {
370  	      SCIPdebugMsg(scip, "over-inequality  %d: %g*x <= %g*y + %g\n", i, overineqs[3*i], overineqs[3*i+1], overineqs[3*i+2]);
371  	      ineqs[3*nineqs] = overineqs[3*i];
372  	      ineqs[3*nineqs+1] = overineqs[3*i+1];
373  	      ineqs[3*nineqs+2] = overineqs[3*i+2];
374  	      ++nineqs;
375  	   }
376  	   for( i = 0; i < nunderineqs; ++i )
377  	   {
378  	      SCIPdebugMsg(scip, "under-inequality %d: %g*x <= %g*y + %g 0\n", i, underineqs[3*i], underineqs[3*i+1], underineqs[3*i+2]);
379  	      ineqs[3*nineqs] = underineqs[3*i];
380  	      ineqs[3*nineqs+1] = underineqs[3*i+1];
381  	      ineqs[3*nineqs+2] = underineqs[3*i+2];
382  	      ++nineqs;
383  	   }
384  	   assert(nineqs == noverineqs + nunderineqs);
385  	
386  	   /* collect children */
387  	   child1 = SCIPexprGetChildren(expr)[0];
388  	   child2 = SCIPexprGetChildren(expr)[1];
389  	   assert(child1 != NULL && child2 != NULL);
390  	   assert(child1 != child2);
391  	
392  	   /* collect bounds of children */
393  	   if( !levelset )
394  	   {
395  	      /* if called from inteval, then use activity */
396  	      boundsx = SCIPexprGetActivity(child1);
397  	      boundsy = SCIPexprGetActivity(child2);
398  	   }
399  	   else
400  	   {
401  	      /* if called from reverseprop, then use bounds */
402  	      boundsx = SCIPgetExprBoundsNonlinear(scip, child1);
403  	      boundsy = SCIPgetExprBoundsNonlinear(scip, child2);
404  	
405  	      /* if children bounds are empty, then returning with *npoints==0 is the way to go */
406  	      if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, boundsx)
407  	          || SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, boundsy) )
408  	         return;
409  	   }
410  	   lbx = boundsx.inf;
411  	   ubx = boundsx.sup;
412  	   lby = boundsy.inf;
413  	   uby = boundsy.sup;
414  	   SCIPdebugMsg(scip, "x = [%g,%g], y=[%g,%g]\n", lbx, ubx, lby, uby);
415  	
416  	   /* corner points that satisfy all inequalities */
417  	   for( i = 0; i < 4; ++i )
418  	   {
419  	      SCIP_Real cx = i < 2 ? lbx : ubx;
420  	      SCIP_Real cy = (i % 2) == 0 ? lby : uby;
421  	
422  	      SCIPdebugMsg(scip, "corner point (%g,%g) feasible? %u\n", cx, cy, isPointFeasible(scip, cx, cy, lbx, ubx, lby, uby, ineqs, nineqs));
423  	
424  	      if( isPointFeasible(scip, cx, cy, lbx, ubx, lby, uby, ineqs, nineqs) )
425  	      {
426  	         xs[*npoints] = cx;
427  	         ys[*npoints] = cy;
428  	         ++(*npoints);
429  	      }
430  	   }
431  	
432  	   /* intersection point of inequalities with [lbx,ubx] x [lby,uby] and extremum of xy on each inequality */
433  	   for( i = 0; i < nineqs; ++i )
434  	   {
435  	      SCIP_Real coefx = ineqs[3*i];
436  	      SCIP_Real coefy = ineqs[3*i+1];
437  	      SCIP_Real constant = ineqs[3*i+2];
438  	      SCIP_Real px[5] = {lbx, ubx, (coefy*lby + constant)/coefx, (coefy*uby + constant)/coefx, 0.0};
439  	      SCIP_Real py[5] = {(coefx*lbx - constant)/coefy, (coefx*ubx - constant)/coefy, lby, uby, 0.0};
440  	      int j;
441  	
442  	      /* the last entry corresponds to the extremum of xy on the line */
443  	      py[4] = (-constant) / (2.0 * coefy);
444  	      px[4] = constant / (2.0 * coefx);
445  	
446  	      for( j = 0; j < 5; ++j )
447  	      {
448  	         SCIPdebugMsg(scip, "intersection point (%g,%g) feasible? %u\n", px[j], py[j], isPointFeasible(scip, px[j], py[j], lbx, ubx, lby, uby, ineqs, nineqs));
449  	         if( isPointFeasible(scip, px[j], py[j], lbx, ubx, lby, uby, ineqs, nineqs) )
450  	         {
451  	            xs[*npoints] = px[j];
452  	            ys[*npoints] = py[j];
453  	            ++(*npoints);
454  	         }
455  	      }
456  	   }
457  	
458  	   /* intersection point between two inequalities */
459  	   for( i = 0; i < nineqs - 1; ++i )
460  	   {
461  	      SCIP_Real coefx1 = ineqs[3*i];
462  	      SCIP_Real coefy1 = ineqs[3*i+1];
463  	      SCIP_Real constant1 = ineqs[3*i+2];
464  	      int j;
465  	
466  	      for( j = i + 1; j < nineqs; ++j )
467  	      {
468  	         SCIP_Real coefx2 = ineqs[3*j];
469  	         SCIP_Real coefy2 = ineqs[3*j+1];
470  	         SCIP_Real constant2 = ineqs[3*j+2];
471  	         SCIP_Real px;
472  	         SCIP_Real py;
473  	
474  	         /* no intersection point -> skip */
475  	         if( SCIPisZero(scip, coefx2*coefy1 - coefx1 * coefy2) )
476  	            continue;
477  	
478  	         py = (constant2 * coefx1 - constant1 * coefx2)/ (coefx2 * coefy1 - coefx1 * coefy2);
479  	         px = (coefy1 * py + constant1) / coefx1;
480  	         assert(SCIPisRelEQ(scip, px, (coefy2 * py + constant2) / coefx2));
481  	
482  	         if( isPointFeasible(scip, px, py, lbx, ubx, lby, uby, ineqs, nineqs) )
483  	         {
484  	            xs[*npoints] = px;
485  	            ys[*npoints] = py;
486  	            ++(*npoints);
487  	         }
488  	      }
489  	   }
490  	
491  	   assert(*npoints <= 22);
492  	
493  	   /* consider the intersection of the level set with
494  	    *
495  	    * 1. the boundary of the box
496  	    * 2. the linear inequalities
497  	    *
498  	    * this adds at most for 4 (level set curves) * 4 (inequalities) * 2 (intersection points) for all linear
499  	    * inequalities and 4 (level set curves) * 2 (intersection points) with the boundary of the box
500  	    */
501  	   if( !levelset )
502  	      return;
503  	
504  	   /* compute intersection of level sets with the boundary */
505  	   for( i = 0; i < 2; ++i )
506  	   {
507  	      SCIP_Real vals[4] = {lbx, ubx, lby, uby};
508  	      SCIP_Real val;
509  	      int k;
510  	
511  	      /* fix auxiliary variable to its lower or upper bound and consider the coefficient of the product */
512  	      val = (i == 0) ? exprbounds.inf : exprbounds.sup;
513  	      val /= SCIPgetCoefExprProduct(expr);
514  	
515  	      for( k = 0; k < 4; ++k )
516  	      {
517  	         if( !SCIPisZero(scip, vals[k]) )
518  	         {
519  	            SCIP_Real res = val / vals[k];
520  	
521  	            assert(SCIPisRelGE(scip, SCIPgetCoefExprProduct(expr)*res*vals[k], exprbounds.inf));
522  	            assert(SCIPisRelLE(scip, SCIPgetCoefExprProduct(expr)*res*vals[k], exprbounds.sup));
523  	
524  	            /* fix x to lbx or ubx */
525  	            if( k < 2 && isPointFeasible(scip, vals[k], res, lbx, ubx, lby, uby, ineqs, nineqs) )
526  	            {
527  	               xs[*npoints] = vals[k];
528  	               ys[*npoints] = res;
529  	               ++(*npoints);
530  	            }
531  	            /* fix y to lby or uby */
532  	            else if( k >= 2 &&  isPointFeasible(scip, res, vals[k], lbx, ubx, lby, uby, ineqs, nineqs) )
533  	            {
534  	               xs[*npoints] = res;
535  	               ys[*npoints] = vals[k];
536  	               ++(*npoints);
537  	            }
538  	         }
539  	      }
540  	   }
541  	
542  	   /* compute intersection points of level sets with the linear inequalities */
543  	   for( i = 0; i < nineqs; ++i )
544  	   {
545  	      SCIP_INTERVAL result;
546  	      SCIP_Real coefx = ineqs[3*i];
547  	      SCIP_Real coefy = ineqs[3*i+1];
548  	      SCIP_Real constant = ineqs[3*i+2];
549  	      SCIP_INTERVAL sqrcoef;
550  	      SCIP_INTERVAL lincoef;
551  	      SCIP_Real px;
552  	      SCIP_Real py;
553  	      int k;
554  	
555  	      /* solve system of coefx x = coefy y + constant and X = xy which is the same as computing the solutions of
556  	       *
557  	       *    (coefy / coefx) y^2 + (constant / coefx) y = inf(X) or sup(X)
558  	       */
559  	      SCIPintervalSet(&sqrcoef, coefy / coefx);
560  	      SCIPintervalSet(&lincoef, constant / coefx);
561  	
562  	      for( k = 0; k < 2; ++k )
563  	      {
564  	         SCIP_INTERVAL rhs;
565  	         SCIP_INTERVAL ybnds;
566  	
567  	         /* set right-hand side */
568  	         if( k == 0 )
569  	            SCIPintervalSet(&rhs, exprbounds.inf);
570  	         else
571  	            SCIPintervalSet(&rhs, exprbounds.sup);
572  	
573  	         SCIPintervalSetBounds(&ybnds, lby, uby);
574  	         SCIPintervalSolveUnivariateQuadExpression(SCIP_INTERVAL_INFINITY, &result, sqrcoef, lincoef, rhs, ybnds);
575  	
576  	         /* interval is empty -> no solution available */
577  	         if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, result) )
578  	            continue;
579  	
580  	         /* compute and check point */
581  	         py = SCIPintervalGetInf(result);
582  	         px = (coefy * py + constant) / coefx;
583  	
584  	         if( isPointFeasible(scip, px, py, lbx, ubx, lby, uby, ineqs, nineqs) )
585  	         {
586  	            xs[*npoints] = px;
587  	            ys[*npoints] = py;
588  	            ++(*npoints);
589  	         }
590  	
591  	         /* check for a second solution */
592  	         if( SCIPintervalGetInf(result) != SCIPintervalGetSup(result) ) /*lint !e777*/
593  	         {
594  	            py = SCIPintervalGetSup(result);
595  	            px = (coefy * py + constant) / coefx;
596  	
597  	            if( isPointFeasible(scip, px, py, lbx, ubx, lby, uby, ineqs, nineqs) )
598  	            {
599  	               xs[*npoints] = px;
600  	               ys[*npoints] = py;
601  	               ++(*npoints);
602  	            }
603  	         }
604  	      }
605  	   }
606  	
607  	   assert(*npoints <= 62);
608  	}
609  	
610  	/** computes interval for a bilinear term when using at least one inequality */
611  	static
612  	SCIP_INTERVAL intevalBilinear(
613  	   SCIP*                 scip,               /**< SCIP data structure */
614  	   SCIP_EXPR*            expr,               /**< product expression */
615  	   SCIP_Real*            underineqs,         /**< inequalities for underestimation */
616  	   int                   nunderineqs,        /**< total number of inequalities for underestimation */
617  	   SCIP_Real*            overineqs,          /**< inequalities for overestimation */
618  	   int                   noverineqs          /**< total number of inequalities for overestimation */
619  	   )
620  	{
621  	   SCIP_INTERVAL interval = {0., 0.};
622  	   SCIP_Real xs[22];
623  	   SCIP_Real ys[22];
624  	   SCIP_Real inf;
625  	   SCIP_Real sup;
626  	   int npoints;
627  	   int i;
628  	
629  	   assert(scip != NULL);
630  	   assert(expr != NULL);
631  	   assert(SCIPexprGetNChildren(expr) == 2);
632  	   assert(noverineqs + nunderineqs <= 4);
633  	
634  	   /* no inequalities available -> skip computation */
635  	   if( noverineqs == 0 && nunderineqs == 0 )
636  	   {
637  	      SCIPintervalSetEntire(SCIP_INTERVAL_INFINITY, &interval);
638  	      return interval;
639  	   }
640  	
641  	   /* x or y has empty interval -> empty */
642  	   if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, SCIPexprGetActivity(SCIPexprGetChildren(expr)[0])) ||
643  	       SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, SCIPexprGetActivity(SCIPexprGetChildren(expr)[1])) )
644  	   {
645  	      SCIPintervalSetEmpty(&interval);
646  	      return interval;
647  	   }
648  	
649  	   /* compute all feasible points (since we use levelset == FALSE, the value of interval doesn't matter) */
650  	   getFeasiblePointsBilinear(scip, NULL, expr, interval, underineqs, nunderineqs, overineqs,
651  	         noverineqs, FALSE, xs, ys, &npoints);
652  	
653  	   /* no feasible point left -> return an empty interval */
654  	   if( npoints == 0 )
655  	   {
656  	      SCIPintervalSetEmpty(&interval);
657  	      return interval;
658  	   }
659  	
660  	   /* compute the minimum and maximum over all computed points */
661  	   inf = xs[0] * ys[0];
662  	   sup = inf;
663  	   SCIPdebugMsg(scip, "point 0: (%g,%g) -> inf = sup = %g\n", xs[0], ys[0], inf);
664  	   for( i = 1; i < npoints; ++i )
665  	   {
666  	      inf = MIN(inf, xs[i] * ys[i]);
667  	      sup = MAX(sup, xs[i] * ys[i]);
668  	      SCIPdebugMsg(scip, "point %d: (%g,%g) -> inf = %g, sup = %g\n", i, xs[i], ys[i], inf, sup);
669  	   }
670  	   assert(inf <= sup);
671  	
672  	   /* adjust infinite values */
673  	   inf = MAX(inf, -SCIP_INTERVAL_INFINITY);
674  	   sup = MIN(sup, SCIP_INTERVAL_INFINITY);
675  	
676  	   /* multiply resulting interval with coefficient of the product expression */
677  	   SCIPintervalSetBounds(&interval, inf, sup);
678  	   if( SCIPgetCoefExprProduct(expr) != 1.0 )
679  	      SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &interval, interval, SCIPgetCoefExprProduct(expr));
680  	
681  	   return interval;
682  	}
683  	
684  	/** uses inequalities for bilinear terms to get stronger bounds during reverse propagation */
685  	static
686  	void reversePropBilinear(
687  	   SCIP*                 scip,               /**< SCIP data structure */
688  	   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler */
689  	   SCIP_EXPR*            expr,               /**< product expression */
690  	   SCIP_INTERVAL         exprbounds,         /**< bounds on product expression */
691  	   SCIP_Real*            underineqs,         /**< inequalities for underestimation */
692  	   int                   nunderineqs,        /**< total number of inequalities for underestimation */
693  	   SCIP_Real*            overineqs,          /**< inequalities for overestimation */
694  	   int                   noverineqs,         /**< total number of inequalities for overestimation */
695  	   SCIP_INTERVAL*        intervalx,          /**< buffer to store the new interval for x */
696  	   SCIP_INTERVAL*        intervaly           /**< buffer to store the new interval for y */
697  	   )
698  	{
699  	   SCIP_Real xs[62];
700  	   SCIP_Real ys[62];
701  	   SCIP_Real exprinf;
702  	   SCIP_Real exprsup;
703  	   SCIP_Bool first = TRUE;
704  	   int npoints;
705  	   int i;
706  	
707  	   assert(scip != NULL);
708  	   assert(conshdlr != NULL);
709  	   assert(expr != NULL);
710  	   assert(intervalx != NULL);
711  	   assert(intervaly != NULL);
712  	   assert(SCIPexprGetNChildren(expr) == 2);
713  	
714  	   assert(noverineqs + nunderineqs > 0);
715  	
716  	   /* set intervals to be empty */
717  	   SCIPintervalSetEmpty(intervalx);
718  	   SCIPintervalSetEmpty(intervaly);
719  	
720  	   /* compute feasible points */
721  	   getFeasiblePointsBilinear(scip, conshdlr, expr, exprbounds, underineqs, nunderineqs, overineqs,
722  	         noverineqs, TRUE, xs, ys, &npoints);
723  	
724  	   /* no feasible points left -> problem is infeasible */
725  	   if( npoints == 0 )
726  	      return;
727  	
728  	   /* get bounds of the product expression */
729  	   exprinf = exprbounds.inf;
730  	   exprsup = exprbounds.sup;
731  	
732  	   /* update intervals with the computed points */
733  	   for( i = 0; i < npoints; ++i )
734  	   {
735  	      SCIP_Real val = SCIPgetCoefExprProduct(expr) * xs[i] * ys[i];
736  	
737  	#ifndef NDEBUG
738  	      {
739  	         SCIP_Real lbx = SCIPexprGetActivity(SCIPexprGetChildren(expr)[0]).inf;
740  	         SCIP_Real ubx = SCIPexprGetActivity(SCIPexprGetChildren(expr)[0]).sup;
741  	         SCIP_Real lby = SCIPexprGetActivity(SCIPexprGetChildren(expr)[1]).inf;
742  	         SCIP_Real uby = SCIPexprGetActivity(SCIPexprGetChildren(expr)[1]).sup;
743  	
744  	         assert(nunderineqs == 0 || isPointFeasible(scip, xs[i], ys[i], lbx, ubx, lby, uby, underineqs, nunderineqs));
745  	         assert(noverineqs == 0 || isPointFeasible(scip, xs[i], ys[i], lbx, ubx, lby, uby, overineqs, noverineqs));
746  	      }
747  	#endif
748  	
749  	      /* only accept points for which the value of x*y is in the interval of the product expression
750  	       *
751  	       * NOTE: in order to consider all relevant points, we are a bit conservative here and relax the interval of
752  	       *      the expression by SCIPfeastol()
753  	       */
754  	      if( SCIPisRelGE(scip, val, exprinf - SCIPfeastol(scip)) && SCIPisRelLE(scip, val, exprsup + SCIPfeastol(scip)) )
755  	      {
756  	         if( first )
757  	         {
758  	            SCIPintervalSet(intervalx, xs[i]);
759  	            SCIPintervalSet(intervaly, ys[i]);
760  	            first = FALSE;
761  	         }
762  	         else
763  	         {
764  	            (*intervalx).inf = MIN((*intervalx).inf, xs[i]);
765  	            (*intervalx).sup = MAX((*intervalx).sup, xs[i]);
766  	            (*intervaly).inf = MIN((*intervaly).inf, ys[i]);
767  	            (*intervaly).sup = MAX((*intervaly).sup, ys[i]);
768  	         }
769  	
770  	         SCIPdebugMsg(scip, "consider points (%g,%g)=%g for reverse propagation\n", xs[i], ys[i], val);
771  	      }
772  	   }
773  	}
774  	
775  	/** helper function to compute the convex envelope of a bilinear term when two linear inequalities are given; we
776  	 *  use the same notation and formulas as in Locatelli 2016
777  	 */
778  	static
779  	void computeBilinEnvelope2(
780  	   SCIP*                 scip,               /**< SCIP data structure */
781  	   SCIP_Real             x,                  /**< reference point for x */
782  	   SCIP_Real             y,                  /**< reference point for y */
783  	   SCIP_Real             mi,                 /**< coefficient of x in the first linear inequality */
784  	   SCIP_Real             qi,                 /**< constant in the first linear inequality */
785  	   SCIP_Real             mj,                 /**< coefficient of x in the second linear inequality */
786  	   SCIP_Real             qj,                 /**< constant in the second linear inequality */
787  	   SCIP_Real* RESTRICT   xi,                 /**< buffer to store x coordinate of the first point */
788  	   SCIP_Real* RESTRICT   yi,                 /**< buffer to store y coordinate of the first point */
789  	   SCIP_Real* RESTRICT   xj,                 /**< buffer to store x coordinate of the second point */
790  	   SCIP_Real* RESTRICT   yj,                 /**< buffer to store y coordinate of the second point */
791  	   SCIP_Real* RESTRICT   xcoef,              /**< buffer to store the x coefficient of the envelope */
792  	   SCIP_Real* RESTRICT   ycoef,              /**< buffer to store the y coefficient of the envelope */
793  	   SCIP_Real* RESTRICT   constant            /**< buffer to store the constant of the envelope */
794  	   )
795  	{
796  	   SCIP_Real QUAD(xiq);
797  	   SCIP_Real QUAD(yiq);
798  	   SCIP_Real QUAD(xjq);
799  	   SCIP_Real QUAD(yjq);
800  	   SCIP_Real QUAD(xcoefq);
801  	   SCIP_Real QUAD(ycoefq);
802  	   SCIP_Real QUAD(constantq);
803  	   SCIP_Real QUAD(tmpq);
804  	
805  	   assert(xi != NULL);
806  	   assert(yi != NULL);
807  	   assert(xj != NULL);
808  	   assert(yj != NULL);
809  	   assert(xcoef != NULL);
810  	   assert(ycoef != NULL);
811  	   assert(constant != NULL);
812  	
813  	   if( SCIPisEQ(scip, mi, mj) )
814  	   {
815  	      /* xi = (x + mi * y - qi) / (2.0*mi) */
816  	      SCIPquadprecProdDD(xiq, mi, y);
817  	      SCIPquadprecSumQD(xiq, xiq, x);
818  	      SCIPquadprecSumQD(xiq, xiq, -qi);
819  	      SCIPquadprecDivQD(xiq, xiq, 2.0 * mi);
820  	      assert(EPSEQ((x + mi * y - qi) / (2.0*mi), QUAD_TO_DBL(xiq), 1e-3));
821  	
822  	      /* yi = mi*(*xi) + qi */
823  	      SCIPquadprecProdQD(yiq, xiq, mi);
824  	      SCIPquadprecSumQD(yiq, yiq, qi);
825  	      assert(EPSEQ(mi*QUAD_TO_DBL(xiq) + qi, QUAD_TO_DBL(yiq), 1e-3));
826  	
827  	      /* xj = (*xi) + (qi - qj)/ (2.0*mi) */
828  	      SCIPquadprecSumDD(xjq, qi, -qj);
829  	      SCIPquadprecDivQD(xjq, xjq, 2.0 * mi);
830  	      SCIPquadprecSumQQ(xjq, xjq, xiq);
831  	      assert(EPSEQ(QUAD_TO_DBL(xiq) + (qi - qj)/ (2.0*mi), QUAD_TO_DBL(xjq), 1e-3));
832  	
833  	      /* yj = mj * (*xj) + qj */
834  	      SCIPquadprecProdQD(yjq, xjq, mj);
835  	      SCIPquadprecSumQD(yjq, yjq, qj);
836  	      assert(EPSEQ(mj * QUAD_TO_DBL(xjq) + qj, QUAD_TO_DBL(yjq), 1e-3));
837  	
838  	      /* ycoef = (*xi) + (qi - qj) / (4.0*mi) note that this is wrong in Locatelli 2016 */
839  	      SCIPquadprecSumDD(ycoefq, qi, -qj);
840  	      SCIPquadprecDivQD(ycoefq, ycoefq, 4.0 * mi);
841  	      SCIPquadprecSumQQ(ycoefq, ycoefq, xiq);
842  	      assert(EPSEQ(QUAD_TO_DBL(xiq) + (qi - qj) / (4.0*mi), QUAD_TO_DBL(ycoefq), 1e-3));
843  	
844  	      /* xcoef = 2.0*mi*(*xi) - mi * (*ycoef) + qi */
845  	      SCIPquadprecProdQD(xcoefq, xiq, 2.0 * mi);
846  	      SCIPquadprecProdQD(tmpq, ycoefq, -mi);
847  	      SCIPquadprecSumQQ(xcoefq, xcoefq, tmpq);
848  	      SCIPquadprecSumQD(xcoefq, xcoefq, qi);
849  	      assert(EPSEQ(2.0*mi*QUAD_TO_DBL(xiq) - mi * QUAD_TO_DBL(ycoefq) + qi, QUAD_TO_DBL(xcoefq), 1e-3));
850  	
851  	      /* constant = -mj*SQR(*xj) - (*ycoef) * qj */
852  	      SCIPquadprecSquareQ(constantq, xjq);
853  	      SCIPquadprecProdQD(constantq, constantq, -mj);
854  	      SCIPquadprecProdQD(tmpq, ycoefq, -qj);
855  	      SCIPquadprecSumQQ(constantq, constantq, tmpq);
856  	      /* assert(EPSEQ(-mj*SQR(QUAD_TO_DBL(xjq)) - QUAD_TO_DBL(ycoefq) * qj, QUAD_TO_DBL(constantq), 1e-3)); */
857  	
858  	      *xi = QUAD_TO_DBL(xiq);
859  	      *yi = QUAD_TO_DBL(yiq);
860  	      *xj = QUAD_TO_DBL(xjq);
861  	      *yj = QUAD_TO_DBL(yjq);
862  	      *ycoef = QUAD_TO_DBL(ycoefq);
863  	      *xcoef = QUAD_TO_DBL(xcoefq);
864  	      *constant = QUAD_TO_DBL(constantq);
865  	   }
866  	   else if( mi > 0.0 )
867  	   {
868  	      assert(mj > 0.0);
869  	
870  	      /* xi = (y + SQRT(mi*mj)*x - qi) / (REALABS(mi) + SQRT(mi*mj)) */
871  	      SCIPquadprecProdDD(xiq, mi, mj);
872  	      SCIPquadprecSqrtQ(xiq, xiq);
873  	      SCIPquadprecProdQD(xiq, xiq, x);
874  	      SCIPquadprecSumQD(xiq, xiq, y);
875  	      SCIPquadprecSumQD(xiq, xiq, -qi); /* (y + SQRT(mi*mj)*x - qi) */
876  	      SCIPquadprecProdDD(tmpq, mi, mj);
877  	      SCIPquadprecSqrtQ(tmpq, tmpq);
878  	      SCIPquadprecSumQD(tmpq, tmpq, REALABS(mi)); /* REALABS(mi) + SQRT(mi*mj) */
879  	      SCIPquadprecDivQQ(xiq, xiq, tmpq);
880  	      assert(EPSEQ((y + SQRT(mi*mj)*x - qi) / (REALABS(mi) + SQRT(mi*mj)), QUAD_TO_DBL(xiq), 1e-3));
881  	
882  	      /* yi = mi*(*xi) + qi */
883  	      SCIPquadprecProdQD(yiq, xiq, mi);
884  	      SCIPquadprecSumQD(yiq, yiq, qi);
885  	      assert(EPSEQ(mi*(QUAD_TO_DBL(xiq)) + qi, QUAD_TO_DBL(yiq), 1e-3));
886  	
887  	      /* xj = (y + SQRT(mi*mj)*x - qj) / (REALABS(mj) + SQRT(mi*mj)) */
888  	      SCIPquadprecProdDD(xjq, mi, mj);
889  	      SCIPquadprecSqrtQ(xjq, xjq);
890  	      SCIPquadprecProdQD(xjq, xjq, x);
891  	      SCIPquadprecSumQD(xjq, xjq, y);
892  	      SCIPquadprecSumQD(xjq, xjq, -qj); /* (y + SQRT(mi*mj)*x - qj) */
893  	      SCIPquadprecProdDD(tmpq, mi, mj);
894  	      SCIPquadprecSqrtQ(tmpq, tmpq);
895  	      SCIPquadprecSumQD(tmpq, tmpq, REALABS(mj)); /* REALABS(mj) + SQRT(mi*mj) */
896  	      SCIPquadprecDivQQ(xjq, xjq, tmpq);
897  	      assert(EPSEQ((y + SQRT(mi*mj)*x - qj) / (REALABS(mj) + SQRT(mi*mj)), QUAD_TO_DBL(xjq), 1e-3));
898  	
899  	      /* yj = mj*(*xj) + qj */
900  	      SCIPquadprecProdQD(yjq, xjq, mj);
901  	      SCIPquadprecSumQD(yjq, yjq, qj);
902  	      assert(EPSEQ(mj*QUAD_TO_DBL(xjq) + qj, QUAD_TO_DBL(yjq), 1e-3));
903  	
904  	      /* ycoef = (2.0*mj*(*xj) + qj - 2.0*mi*(*xi) - qi) / (mj - mi) */
905  	      SCIPquadprecProdQD(ycoefq, xjq, 2.0 * mj);
906  	      SCIPquadprecSumQD(ycoefq, ycoefq, qj);
907  	      SCIPquadprecProdQD(tmpq, xiq, -2.0 * mi);
908  	      SCIPquadprecSumQQ(ycoefq, ycoefq, tmpq);
909  	      SCIPquadprecSumQD(ycoefq, ycoefq, -qi);
910  	      SCIPquadprecSumDD(tmpq, mj, -mi);
911  	      SCIPquadprecDivQQ(ycoefq, ycoefq, tmpq);
912  	      assert(EPSEQ((2.0*mj*QUAD_TO_DBL(xjq) + qj - 2.0*mi*QUAD_TO_DBL(xiq) - qi) / (mj - mi), QUAD_TO_DBL(ycoefq), 1e-3));
913  	
914  	      /* xcoef = 2.0*mj*(*xj) + qj - mj*(*ycoef) */
915  	      SCIPquadprecProdQD(xcoefq, xjq, 2.0 * mj);
916  	      SCIPquadprecSumQD(xcoefq, xcoefq, qj);
917  	      SCIPquadprecProdQD(tmpq, ycoefq, -mj);
918  	      SCIPquadprecSumQQ(xcoefq, xcoefq, tmpq);
919  	      assert(EPSEQ(2.0*mj*QUAD_TO_DBL(xjq) + qj - mj*QUAD_TO_DBL(ycoefq), QUAD_TO_DBL(xcoefq), 1e-3));
920  	
921  	      /* constant = -mj*SQR(*xj) - (*ycoef) * qj */
922  	      SCIPquadprecSquareQ(constantq, xjq);
923  	      SCIPquadprecProdQD(constantq, constantq, -mj);
924  	      SCIPquadprecProdQD(tmpq, ycoefq, -qj);
925  	      SCIPquadprecSumQQ(constantq, constantq, tmpq);
926  	      /* assert(EPSEQ(-mj*SQR(QUAD_TO_DBL(xjq)) - QUAD_TO_DBL(ycoefq) * qj, QUAD_TO_DBL(constantq), 1e-3)); */
927  	
928  	      *xi = QUAD_TO_DBL(xiq);
929  	      *yi = QUAD_TO_DBL(yiq);
930  	      *xj = QUAD_TO_DBL(xjq);
931  	      *yj = QUAD_TO_DBL(yjq);
932  	      *ycoef = QUAD_TO_DBL(ycoefq);
933  	      *xcoef = QUAD_TO_DBL(xcoefq);
934  	      *constant = QUAD_TO_DBL(constantq);
935  	   }
936  	   else
937  	   {
938  	      assert(mi < 0.0 && mj < 0.0);
939  	
940  	      /* apply variable transformation x = -x in case for overestimation */
941  	      computeBilinEnvelope2(scip, -x, y, -mi, qi, -mj, qj, xi, yi, xj, yj, xcoef, ycoef, constant);
942  	
943  	      /* revert transformation; multiply cut by -1 and change -x by x */
944  	      *xi = -(*xi);
945  	      *xj = -(*xj);
946  	      *ycoef = -(*ycoef);
947  	      *constant = -(*constant);
948  	   }
949  	}
950  	
951  	/** output method of statistics table to output file stream 'file' */
952  	static
953  	SCIP_DECL_TABLEOUTPUT(tableOutputBilinear)
954  	{ /*lint --e{715}*/
955  	   SCIP_NLHDLR* nlhdlr;
956  	   SCIP_NLHDLRDATA* nlhdlrdata;
957  	   SCIP_CONSHDLR* conshdlr;
958  	   SCIP_HASHMAP* hashmap;
959  	   SCIP_EXPRITER* it;
960  	   int resfound = 0;
961  	   int restotal = 0;
962  	   int c;
963  	
964  	   conshdlr = SCIPfindConshdlr(scip, "nonlinear");
965  	   assert(conshdlr != NULL);
966  	   nlhdlr = SCIPfindNlhdlrNonlinear(conshdlr, NLHDLR_NAME);
967  	   assert(nlhdlr != NULL);
968  	   nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
969  	   assert(nlhdlrdata != NULL);
970  	
971  	   /* allocate memory */
972  	   SCIP_CALL( SCIPhashmapCreate(&hashmap, SCIPblkmem(scip), nlhdlrdata->nexprs) );
973  	   SCIP_CALL( SCIPcreateExpriter(scip, &it) );
974  	
975  	   for( c = 0; c < nlhdlrdata->nexprs; ++c )
976  	   {
977  	      assert(!SCIPhashmapExists(hashmap, nlhdlrdata->exprs[c]));
978  	      SCIP_CALL( SCIPhashmapInsertInt(hashmap, nlhdlrdata->exprs[c], 0) );
979  	   }
980  	
981  	   /* count in how many constraints each expression is contained */
982  	   for( c = 0; c < SCIPconshdlrGetNConss(conshdlr); ++c )
983  	   {
984  	      SCIP_CONS* cons = SCIPconshdlrGetConss(conshdlr)[c];
985  	      SCIP_EXPR* expr;
986  	
987  	      SCIP_CALL( SCIPexpriterInit(it, SCIPgetExprNonlinear(cons), SCIP_EXPRITER_DFS, FALSE) );
988  	
989  	      for( expr = SCIPexpriterGetCurrent(it); !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) ) /*lint !e441*//*lint !e440*/
990  	      {
991  	         if( SCIPhashmapExists(hashmap, expr) )
992  	         {
993  	            int nuses = SCIPhashmapGetImageInt(hashmap, expr);
994  	            SCIP_CALL( SCIPhashmapSetImageInt(hashmap, expr, nuses + 1) );
995  	         }
996  	      }
997  	   }
998  	
999  	   /* compute success ratio */
1000 	   for( c = 0; c < nlhdlrdata->nexprs; ++c )
1001 	   {
1002 	      SCIP_NLHDLREXPRDATA* nlhdlrexprdata;
1003 	      int nuses;
1004 	
1005 	      nuses = SCIPhashmapGetImageInt(hashmap, nlhdlrdata->exprs[c]);
1006 	      assert(nuses > 0);
1007 	
1008 	      nlhdlrexprdata = SCIPgetNlhdlrExprDataNonlinear(nlhdlr, nlhdlrdata->exprs[c]);
1009 	      assert(nlhdlrexprdata != NULL);
1010 	
1011 	      if( nlhdlrexprdata->nunderineqs > 0 || nlhdlrexprdata->noverineqs > 0 )
1012 	         resfound += nuses;
1013 	      restotal += nuses;
1014 	   }
1015 	
1016 	   /* print statistics */
1017 	   SCIPinfoMessage(scip, file, "Bilinear Nlhdlr    : %10s %10s\n", "#found", "#total");
1018 	   SCIPinfoMessage(scip, file, "  %-17s:", "");
1019 	   SCIPinfoMessage(scip, file, " %10d", resfound);
1020 	   SCIPinfoMessage(scip, file, " %10d", restotal);
1021 	   SCIPinfoMessage(scip, file, "\n");
1022 	
1023 	   /* free memory */
1024 	   SCIPfreeExpriter(&it);
1025 	   SCIPhashmapFree(&hashmap);
1026 	
1027 	   return SCIP_OKAY;
1028 	}
1029 	
1030 	
1031 	/*
1032 	 * Callback methods of nonlinear handler
1033 	 */
1034 	
1035 	/** nonlinear handler copy callback */
1036 	static
1037 	SCIP_DECL_NLHDLRCOPYHDLR(nlhdlrCopyhdlrBilinear)
1038 	{ /*lint --e{715}*/
1039 	   assert(targetscip != NULL);
1040 	   assert(sourcenlhdlr != NULL);
1041 	   assert(strcmp(SCIPnlhdlrGetName(sourcenlhdlr), NLHDLR_NAME) == 0);
1042 	
1043 	   SCIP_CALL( SCIPincludeNlhdlrBilinear(targetscip) );
1044 	
1045 	   return SCIP_OKAY;
1046 	}
1047 	
1048 	/** callback to free data of handler */
1049 	static
1050 	SCIP_DECL_NLHDLRFREEHDLRDATA(nlhdlrFreehdlrdataBilinear)
1051 	{ /*lint --e{715}*/
1052 	   assert(nlhdlrdata != NULL);
1053 	   assert((*nlhdlrdata)->nexprs == 0);
1054 	
1055 	   if( (*nlhdlrdata)->exprmap != NULL )
1056 	   {
1057 	      assert(SCIPhashmapGetNElements((*nlhdlrdata)->exprmap) == 0);
1058 	      SCIPhashmapFree(&(*nlhdlrdata)->exprmap);
1059 	   }
1060 	
1061 	   SCIPfreeBlockMemoryArrayNull(scip, &(*nlhdlrdata)->exprs, (*nlhdlrdata)->exprsize);
1062 	   SCIPfreeBlockMemory(scip, nlhdlrdata);
1063 	
1064 	   return SCIP_OKAY;
1065 	}
1066 	
1067 	/** callback to free expression specific data */
1068 	static
1069 	SCIP_DECL_NLHDLRFREEEXPRDATA(nlhdlrFreeExprDataBilinear)
1070 	{  /*lint --e{715}*/
1071 	   SCIP_NLHDLRDATA* nlhdlrdata;
1072 	   int pos;
1073 	
1074 	   assert(expr != NULL);
1075 	
1076 	   nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
1077 	   assert(nlhdlrdata != NULL);
1078 	   assert(nlhdlrdata->nexprs > 0);
1079 	   assert(nlhdlrdata->exprs != NULL);
1080 	   assert(nlhdlrdata->exprmap != NULL);
1081 	   assert(SCIPhashmapExists(nlhdlrdata->exprmap, (void*)expr));
1082 	
1083 	   pos = SCIPhashmapGetImageInt(nlhdlrdata->exprmap, (void*)expr);
1084 	   assert(pos >= 0 && pos < nlhdlrdata->nexprs);
1085 	   assert(nlhdlrdata->exprs[pos] == expr);
1086 	
1087 	   /* move the last expression to the free position */
1088 	   if( nlhdlrdata->nexprs > 0 && pos != nlhdlrdata->nexprs - 1 )
1089 	   {
1090 	      SCIP_EXPR* lastexpr = nlhdlrdata->exprs[nlhdlrdata->nexprs - 1];
1091 	      assert(expr != lastexpr);
1092 	      assert(SCIPhashmapExists(nlhdlrdata->exprmap, (void*)lastexpr));
1093 	
1094 	      nlhdlrdata->exprs[pos] = lastexpr;
1095 	      nlhdlrdata->exprs[nlhdlrdata->nexprs - 1] = NULL;
1096 	      SCIP_CALL( SCIPhashmapSetImageInt(nlhdlrdata->exprmap, (void*)lastexpr, pos) );
1097 	   }
1098 	
1099 	   /* remove expression from the nonlinear handler data */
1100 	   SCIP_CALL( SCIPhashmapRemove(nlhdlrdata->exprmap, (void*)expr) );
1101 	   SCIP_CALL( SCIPreleaseExpr(scip, &expr) );
1102 	   --nlhdlrdata->nexprs;
1103 	
1104 	   /* free nonlinear handler expression data */
1105 	   SCIPfreeBlockMemoryNull(scip, nlhdlrexprdata);
1106 	
1107 	   return SCIP_OKAY;
1108 	}
1109 	
1110 	/** callback to be called in initialization */
1111 	#define nlhdlrInitBilinear NULL
1112 	
1113 	/** callback to be called in deinitialization */
1114 	static
1115 	SCIP_DECL_NLHDLREXIT(nlhdlrExitBilinear)
1116 	{  /*lint --e{715}*/
1117 	   assert(SCIPnlhdlrGetData(nlhdlr) != NULL);
1118 	   assert(SCIPnlhdlrGetData(nlhdlr)->nexprs == 0);
1119 	
1120 	   return SCIP_OKAY;
1121 	}
1122 	
1123 	/** callback to detect structure in expression tree */
1124 	static
1125 	SCIP_DECL_NLHDLRDETECT(nlhdlrDetectBilinear)
1126 	{ /*lint --e{715}*/
1127 	   SCIP_NLHDLRDATA* nlhdlrdata;
1128 	
1129 	   assert(expr != NULL);
1130 	   assert(participating != NULL);
1131 	
1132 	   nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
1133 	   assert(nlhdlrdata);
1134 	
1135 	   /* only during solving will we have the extra inequalities that we rely on so much here */
1136 	   if( SCIPgetStage(scip) < SCIP_STAGE_INITSOLVE )
1137 	      return SCIP_OKAY;
1138 	
1139 	   /* check for product expressions with two children */
1140 	   if( SCIPisExprProduct(scip, expr) && SCIPexprGetNChildren(expr) == 2
1141 	      && (nlhdlrdata->exprmap == NULL || !SCIPhashmapExists(nlhdlrdata->exprmap, (void*)expr)) )
1142 	   {
1143 	      SCIP_EXPR** children;
1144 	      SCIP_Bool valid;
1145 	      int c;
1146 	
1147 	      children = SCIPexprGetChildren(expr);
1148 	      assert(children != NULL);
1149 	
1150 	      /* detection is only successful if both children will have auxiliary variable or are variables
1151 	       * that are not binary variables */
1152 	      valid = TRUE;
1153 	      for( c = 0; c < 2; ++c )
1154 	      {
1155 	         assert(children[c] != NULL);
1156 	         if( SCIPgetExprNAuxvarUsesNonlinear(children[c]) == 0 &&
1157 	            (!SCIPisExprVar(scip, children[c]) || SCIPvarIsBinary(SCIPgetVarExprVar(children[c]))) )
1158 	         {
1159 	            valid = FALSE;
1160 	            break;
1161 	         }
1162 	      }
1163 	
1164 	      if( valid )
1165 	      {
1166 	         /* create expression data for the nonlinear handler */
1167 	         SCIP_CALL( SCIPallocClearBlockMemory(scip, nlhdlrexprdata) );
1168 	         (*nlhdlrexprdata)->lastnodeid = -1;
1169 	
1170 	         /* ensure that there is enough memory to store the detected expression */
1171 	         if( nlhdlrdata->exprsize < nlhdlrdata->nexprs + 1 )
1172 	         {
1173 	            int newsize = SCIPcalcMemGrowSize(scip, nlhdlrdata->nexprs + 1);
1174 	            assert(newsize > nlhdlrdata->exprsize);
1175 	
1176 	            SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &nlhdlrdata->exprs, nlhdlrdata->exprsize, newsize) );
1177 	            nlhdlrdata->exprsize = newsize;
1178 	         }
1179 	
1180 	         /* create expression map, if not done so far */
1181 	         if( nlhdlrdata->exprmap == NULL )
1182 	         {
1183 	            SCIP_CALL( SCIPhashmapCreate(&nlhdlrdata->exprmap, SCIPblkmem(scip), SCIPgetNVars(scip)) );
1184 	         }
1185 	
1186 	#ifndef NDEBUG
1187 	         {
1188 	            int i;
1189 	
1190 	            for( i = 0; i < nlhdlrdata->nexprs; ++i )
1191 	               assert(nlhdlrdata->exprs[i] != expr);
1192 	         }
1193 	#endif
1194 	
1195 	         /* add expression to nlhdlrdata and capture it */
1196 	         nlhdlrdata->exprs[nlhdlrdata->nexprs] = expr;
1197 	         SCIPcaptureExpr(expr);
1198 	         SCIP_CALL( SCIPhashmapInsertInt(nlhdlrdata->exprmap, (void*)expr, nlhdlrdata->nexprs) );
1199 	         ++nlhdlrdata->nexprs;
1200 	
1201 	         /* tell children that we will use their auxvar and use its activity for both estimate and domain propagation */
1202 	         SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, children[0], TRUE, nlhdlrdata->useinteval
1203 	               || nlhdlrdata->usereverseprop, TRUE, TRUE) );
1204 	         SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, children[1], TRUE, nlhdlrdata->useinteval
1205 	               || nlhdlrdata->usereverseprop, TRUE, TRUE) );
1206 	      }
1207 	   }
1208 	
1209 	   if( *nlhdlrexprdata != NULL )
1210 	   {
1211 	      /* we want to join separation and domain propagation, if not disabled by parameter */
1212 	      *participating = SCIP_NLHDLR_METHOD_SEPABOTH;
1213 	      if( nlhdlrdata->useinteval || nlhdlrdata->usereverseprop )
1214 	         *participating |= SCIP_NLHDLR_METHOD_ACTIVITY;
1215 	   }
1216 	
1217 	#ifdef SCIP_DEBUG
1218 	   if( *participating )
1219 	   {
1220 	      SCIPdebugMsg(scip, "detected expr ");
1221 	      SCIPprintExpr(scip, expr, NULL);
1222 	      SCIPinfoMessage(scip, NULL, " participating: %d\n", *participating);
1223 	   }
1224 	#endif
1225 	
1226 	   return SCIP_OKAY;
1227 	}
1228 	
1229 	/** auxiliary evaluation callback of nonlinear handler */
1230 	static
1231 	SCIP_DECL_NLHDLREVALAUX(nlhdlrEvalauxBilinear)
1232 	{ /*lint --e{715}*/
1233 	   SCIP_VAR* var1;
1234 	   SCIP_VAR* var2;
1235 	   SCIP_Real coef;
1236 	
1237 	   assert(SCIPisExprProduct(scip, expr));
1238 	   assert(SCIPexprGetNChildren(expr) == 2);
1239 	
1240 	   var1 = SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(expr)[0]);
1241 	   assert(var1 != NULL);
1242 	   var2 = SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(expr)[1]);
1243 	   assert(var2 != NULL);
1244 	   coef = SCIPgetCoefExprProduct(expr);
1245 	
1246 	   *auxvalue = coef * SCIPgetSolVal(scip, sol, var1) * SCIPgetSolVal(scip, sol, var2);
1247 	
1248 	   return SCIP_OKAY;
1249 	}
1250 	
1251 	/** separation initialization method of a nonlinear handler (called during CONSINITLP) */
1252 	#define nlhdlrInitSepaBilinear NULL
1253 	
1254 	/** separation deinitialization method of a nonlinear handler (called during CONSEXITSOL) */
1255 	#define nlhdlrExitSepaBilinear NULL
1256 	
1257 	/** nonlinear handler separation callback */
1258 	#define nlhdlrEnfoBilinear NULL
1259 	
1260 	/** nonlinear handler under/overestimation callback */
1261 	static
1262 	SCIP_DECL_NLHDLRESTIMATE(nlhdlrEstimateBilinear)
1263 	{ /*lint --e{715}*/
1264 	   SCIP_NLHDLRDATA* nlhdlrdata;
1265 	   SCIP_VAR* x;
1266 	   SCIP_VAR* y;
1267 	   SCIP_VAR* auxvar;
1268 	   SCIP_Real lincoefx = 0.0;
1269 	   SCIP_Real lincoefy = 0.0;
1270 	   SCIP_Real linconstant = 0.0;
1271 	   SCIP_Real refpointx;
1272 	   SCIP_Real refpointy;
1273 	   SCIP_Real violation;
1274 	   SCIP_Longint nodeid;
1275 	   SCIP_Bool mccsuccess = TRUE;
1276 	   SCIP_ROWPREP* rowprep;
1277 	
1278 	   assert(rowpreps != NULL);
1279 	
1280 	   *success = FALSE;
1281 	   *addedbranchscores = FALSE;
1282 	
1283 	   /* check whether an inequality is available */
1284 	   if( nlhdlrexprdata->noverineqs == 0 && nlhdlrexprdata->nunderineqs == 0 )
1285 	      return SCIP_OKAY;
1286 	
1287 	   nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
1288 	   assert(nlhdlrdata != NULL);
1289 	
1290 	   nodeid = SCIPnodeGetNumber(SCIPgetCurrentNode(scip));
1291 	
1292 	   /* update last node */
1293 	   if( nlhdlrexprdata->lastnodeid != nodeid )
1294 	   {
1295 	      nlhdlrexprdata->lastnodeid = nodeid;
1296 	      nlhdlrexprdata->nseparoundslastnode = 0;
1297 	   }
1298 	
1299 	   /* update separation round */
1300 	   ++nlhdlrexprdata->nseparoundslastnode;
1301 	
1302 	   /* check working limits */
1303 	   if( (SCIPgetDepth(scip) == 0 && nlhdlrexprdata->nseparoundslastnode > nlhdlrdata->maxseparoundsroot)
1304 	         || (SCIPgetDepth(scip) > 0 && nlhdlrexprdata->nseparoundslastnode > nlhdlrdata->maxseparounds)
1305 	         || SCIPgetDepth(scip) > nlhdlrdata->maxsepadepth )
1306 	      return SCIP_OKAY;
1307 	
1308 	   /* collect variables */
1309 	   x = SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(expr)[0]);
1310 	   assert(x != NULL);
1311 	   y = SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(expr)[1]);
1312 	   assert(y != NULL);
1313 	   auxvar = SCIPgetExprAuxVarNonlinear(expr);
1314 	   assert(auxvar != NULL);
1315 	
1316 	   /* get and adjust the reference points */
1317 	   refpointx = MIN(MAX(SCIPgetSolVal(scip, sol, x), SCIPvarGetLbLocal(x)),SCIPvarGetUbLocal(x)); /*lint !e666*/
1318 	   refpointy = MIN(MAX(SCIPgetSolVal(scip, sol, y), SCIPvarGetLbLocal(y)),SCIPvarGetUbLocal(y)); /*lint !e666*/
1319 	   assert(SCIPisLE(scip, refpointx, SCIPvarGetUbLocal(x)) && SCIPisGE(scip, refpointx, SCIPvarGetLbLocal(x)));
1320 	   assert(SCIPisLE(scip, refpointy, SCIPvarGetUbLocal(y)) && SCIPisGE(scip, refpointy, SCIPvarGetLbLocal(y)));
1321 	
1322 	   /* use McCormick inequalities to decide whether we want to separate or not */
1323 	   SCIPaddBilinMcCormick(scip, SCIPgetCoefExprProduct(expr), SCIPvarGetLbLocal(x), SCIPvarGetUbLocal(x), refpointx,
1324 	         SCIPvarGetLbLocal(y), SCIPvarGetUbLocal(y), refpointy, overestimate, &lincoefx, &lincoefy, &linconstant,
1325 	         &mccsuccess);
1326 	
1327 	   /* too large values in McCormick inequalities -> skip */
1328 	   if( !mccsuccess )
1329 	      return SCIP_OKAY;
1330 	
1331 	   /* compute violation for the McCormick relaxation */
1332 	   violation = lincoefx * refpointx + lincoefy * refpointy + linconstant - SCIPgetSolVal(scip, sol, auxvar);
1333 	   if( overestimate )
1334 	      violation = -violation;
1335 	
1336 	   /* only use a tighter relaxations if McCormick does not separate the reference point */
1337 	   if( SCIPisFeasLE(scip, violation, 0.0) && useBilinIneqs(scip, x, y, refpointx, refpointy) )
1338 	   {
1339 	      SCIP_Bool useoverestineq = SCIPgetCoefExprProduct(expr) > 0.0 ? overestimate : !overestimate;
1340 	      SCIP_Real mccormickval = lincoefx * refpointx + lincoefy * refpointy + linconstant;
1341 	      SCIP_Real* ineqs;
1342 	      SCIP_Real bestval;
1343 	      int nineqs;
1344 	
1345 	      /* McCormick relaxation is too weak */
1346 	      bestval = mccormickval;
1347 	
1348 	      /* get the inequalities that might lead to a tighter relaxation */
1349 	      if( useoverestineq )
1350 	      {
1351 	         ineqs = nlhdlrexprdata->overineqs;
1352 	         nineqs = nlhdlrexprdata->noverineqs;
1353 	      }
1354 	      else
1355 	      {
1356 	         ineqs = nlhdlrexprdata->underineqs;
1357 	         nineqs = nlhdlrexprdata->nunderineqs;
1358 	      }
1359 	
1360 	      /* use linear inequalities to update relaxation */
1361 	      updateBilinearRelaxation(scip, x, y, SCIPgetCoefExprProduct(expr),
1362 	         overestimate ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT,
1363 	         refpointx, refpointy, ineqs, nineqs, mccormickval,
1364 	         &lincoefx, &lincoefy, &linconstant, &bestval,
1365 	         success);
1366 	
1367 	#ifndef NDEBUG
1368 	      /* check whether cut is really valid */
1369 	      if( *success )
1370 	      {
1371 	         assert(!overestimate || SCIPisLE(scip, bestval, mccormickval));
1372 	         assert(overestimate || SCIPisGE(scip, bestval, mccormickval));
1373 	      }
1374 	#endif
1375 	   }
1376 	
1377 	   if( *success )
1378 	   {
1379 	      SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, overestimate ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT, TRUE) );
1380 	      SCIProwprepAddConstant(rowprep, linconstant);
1381 	      SCIP_CALL( SCIPensureRowprepSize(scip, rowprep, 2) );
1382 	      SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, x, lincoefx) );
1383 	      SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, y, lincoefy) );
1384 	      SCIP_CALL( SCIPsetPtrarrayVal(scip, rowpreps, 0, rowprep) );
1385 	   }
1386 	
1387 	   return SCIP_OKAY;
1388 	}
1389 	
1390 	/** nonlinear handler interval evaluation callback */
1391 	static
1392 	SCIP_DECL_NLHDLRINTEVAL(nlhdlrIntevalBilinear)
1393 	{ /*lint --e{715}*/
1394 	   SCIP_NLHDLRDATA* nlhdlrdata;
1395 	   assert(nlhdlrexprdata != NULL);
1396 	
1397 	   nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
1398 	   assert(nlhdlrdata != NULL);
1399 	
1400 	   if( nlhdlrdata->useinteval && nlhdlrexprdata->nunderineqs + nlhdlrexprdata->noverineqs > 0 )
1401 	   {
1402 	      SCIP_INTERVAL tmp = intevalBilinear(scip, expr, nlhdlrexprdata->underineqs, nlhdlrexprdata->nunderineqs,
1403 	         nlhdlrexprdata->overineqs, nlhdlrexprdata->noverineqs);
1404 	
1405 	      /* intersect intervals if we have learned a tighter interval */
1406 	      if( SCIPisGT(scip, tmp.inf, (*interval).inf) || SCIPisLT(scip, tmp.sup, (*interval).sup) )
1407 	         SCIPintervalIntersect(interval, *interval, tmp);
1408 	   }
1409 	
1410 	   return SCIP_OKAY;
1411 	}
1412 	
1413 	/** nonlinear handler callback for reverse propagation */
1414 	static
1415 	SCIP_DECL_NLHDLRREVERSEPROP(nlhdlrReversepropBilinear)
1416 	{ /*lint --e{715}*/
1417 	   SCIP_NLHDLRDATA* nlhdlrdata;
1418 	
1419 	   assert(nlhdlrexprdata != NULL);
1420 	
1421 	   nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
1422 	   assert(nlhdlrdata != NULL);
1423 	
1424 	   if( nlhdlrdata->usereverseprop && nlhdlrexprdata->nunderineqs + nlhdlrexprdata->noverineqs > 0 )
1425 	   {
1426 	      SCIP_EXPR* childx;
1427 	      SCIP_EXPR* childy;
1428 	      SCIP_INTERVAL intervalx;
1429 	      SCIP_INTERVAL intervaly;
1430 	
1431 	      assert(SCIPexprGetNChildren(expr) == 2);
1432 	      childx = SCIPexprGetChildren(expr)[0];
1433 	      childy = SCIPexprGetChildren(expr)[1];
1434 	      assert(childx != NULL && childy != NULL);
1435 	
1436 	      SCIPintervalSetEntire(SCIP_INTERVAL_INFINITY, &intervalx);
1437 	      SCIPintervalSetEntire(SCIP_INTERVAL_INFINITY, &intervaly);
1438 	
1439 	      /* compute bounds on x and y */
1440 	      reversePropBilinear(scip, conshdlr, expr, bounds, nlhdlrexprdata->underineqs, nlhdlrexprdata->nunderineqs,
1441 	         nlhdlrexprdata->overineqs, nlhdlrexprdata->noverineqs, &intervalx, &intervaly);
1442 	
1443 	      /* tighten bounds of x */
1444 	      SCIPdebugMsg(scip, "try to tighten bounds of x: [%g,%g] -> [%g,%g]\n",
1445 	         SCIPgetExprBoundsNonlinear(scip, childx).inf, SCIPgetExprBoundsNonlinear(scip, childx).sup,
1446 	         intervalx.inf, intervalx.sup);
1447 	
1448 	      SCIP_CALL( SCIPtightenExprIntervalNonlinear(scip, SCIPexprGetChildren(expr)[0], intervalx, infeasible,
1449 	         nreductions) );
1450 	
1451 	      if( !(*infeasible) )
1452 	      {
1453 	         /* tighten bounds of y */
1454 	         SCIPdebugMsg(scip, "try to tighten bounds of y: [%g,%g] -> [%g,%g]\n",
1455 	            SCIPgetExprBoundsNonlinear(scip, childy).inf, SCIPgetExprBoundsNonlinear(scip, childy).sup,
1456 	            intervaly.inf, intervaly.sup);
1457 	         SCIP_CALL( SCIPtightenExprIntervalNonlinear(scip, SCIPexprGetChildren(expr)[1], intervaly,
1458 	            infeasible, nreductions) );
1459 	      }
1460 	   }
1461 	
1462 	   return SCIP_OKAY;
1463 	}
1464 	
1465 	/*
1466 	 * nonlinear handler specific interface methods
1467 	 */
1468 	
1469 	/** includes bilinear nonlinear handler in nonlinear constraint handler */
1470 	SCIP_RETCODE SCIPincludeNlhdlrBilinear(
1471 	   SCIP*                 scip                /**< SCIP data structure */
1472 	   )
1473 	{
1474 	   SCIP_NLHDLRDATA* nlhdlrdata;
1475 	   SCIP_NLHDLR* nlhdlr;
1476 	
1477 	   assert(scip != NULL);
1478 	
1479 	   /**! [SnippetIncludeNlhdlrBilinear] */
1480 	   /* create nonlinear handler specific data */
1481 	   SCIP_CALL( SCIPallocBlockMemory(scip, &nlhdlrdata) );
1482 	   BMSclearMemory(nlhdlrdata);
1483 	
1484 	   SCIP_CALL( SCIPincludeNlhdlrNonlinear(scip, &nlhdlr, NLHDLR_NAME, NLHDLR_DESC, NLHDLR_DETECTPRIORITY,
1485 	      NLHDLR_ENFOPRIORITY, nlhdlrDetectBilinear, nlhdlrEvalauxBilinear, nlhdlrdata) );
1486 	   assert(nlhdlr != NULL);
1487 	
1488 	   SCIPnlhdlrSetCopyHdlr(nlhdlr, nlhdlrCopyhdlrBilinear);
1489 	   SCIPnlhdlrSetFreeHdlrData(nlhdlr, nlhdlrFreehdlrdataBilinear);
1490 	   SCIPnlhdlrSetFreeExprData(nlhdlr, nlhdlrFreeExprDataBilinear);
1491 	   SCIPnlhdlrSetInitExit(nlhdlr, nlhdlrInitBilinear, nlhdlrExitBilinear);
1492 	   SCIPnlhdlrSetSepa(nlhdlr, nlhdlrInitSepaBilinear, nlhdlrEnfoBilinear, nlhdlrEstimateBilinear, nlhdlrExitSepaBilinear);
1493 	   SCIPnlhdlrSetProp(nlhdlr, nlhdlrIntevalBilinear, nlhdlrReversepropBilinear);
1494 	
1495 	   /* parameters */
1496 	   SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/useinteval",
1497 	         "whether to use the interval evaluation callback of the nlhdlr",
1498 	         &nlhdlrdata->useinteval, FALSE, TRUE, NULL, NULL) );
1499 	
1500 	   SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/usereverseprop",
1501 	         "whether to use the reverse propagation callback of the nlhdlr",
1502 	         &nlhdlrdata->usereverseprop, FALSE, TRUE, NULL, NULL) );
1503 	
1504 	   SCIP_CALL( SCIPaddIntParam(scip, "nlhdlr/" NLHDLR_NAME "/maxseparoundsroot",
1505 	         "maximum number of separation rounds in the root node",
1506 	         &nlhdlrdata->maxseparoundsroot, FALSE, 10, 0, INT_MAX, NULL, NULL) );
1507 	
1508 	   SCIP_CALL( SCIPaddIntParam(scip, "nlhdlr/" NLHDLR_NAME "/maxseparounds",
1509 	         "maximum number of separation rounds in a local node",
1510 	         &nlhdlrdata->maxseparounds, FALSE, 1, 0, INT_MAX, NULL, NULL) );
1511 	
1512 	   SCIP_CALL( SCIPaddIntParam(scip, "nlhdlr/" NLHDLR_NAME "/maxsepadepth",
1513 	         "maximum depth to apply separation",
1514 	         &nlhdlrdata->maxsepadepth, FALSE, INT_MAX, 0, INT_MAX, NULL, NULL) );
1515 	
1516 	   /* statistic table */
1517 	   assert(SCIPfindTable(scip, TABLE_NAME_BILINEAR) == NULL);
1518 	   SCIP_CALL( SCIPincludeTable(scip, TABLE_NAME_BILINEAR, TABLE_DESC_BILINEAR, FALSE,
1519 	         NULL, NULL, NULL, NULL, NULL, NULL, tableOutputBilinear,
1520 	         NULL, TABLE_POSITION_BILINEAR, TABLE_EARLIEST_STAGE_BILINEAR) );
1521 	   /**! [SnippetIncludeNlhdlrBilinear] */
1522 	
1523 	   return SCIP_OKAY;
1524 	}
1525 	
1526 	/** returns an array of expressions that have been detected by the bilinear nonlinear handler */
1527 	SCIP_EXPR** SCIPgetExprsBilinear(
1528 	   SCIP_NLHDLR*          nlhdlr              /**< nonlinear handler */
1529 	   )
1530 	{
1531 	   SCIP_NLHDLRDATA* nlhdlrdata;
1532 	
1533 	   assert(nlhdlr != NULL);
1534 	   assert(strcmp(SCIPnlhdlrGetName(nlhdlr), NLHDLR_NAME) == 0);
1535 	
1536 	   nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
1537 	   assert(nlhdlrdata);
1538 	
1539 	   return nlhdlrdata->exprs;
1540 	}
1541 	
1542 	/** returns the total number of expressions that have been detected by the bilinear nonlinear handler */
1543 	int SCIPgetNExprsBilinear(
1544 	   SCIP_NLHDLR*          nlhdlr              /**< nonlinear handler */
1545 	   )
1546 	{
1547 	   SCIP_NLHDLRDATA* nlhdlrdata;
1548 	
1549 	   assert(nlhdlr != NULL);
1550 	   assert(strcmp(SCIPnlhdlrGetName(nlhdlr), NLHDLR_NAME) == 0);
1551 	
1552 	   nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
1553 	   assert(nlhdlrdata);
1554 	
1555 	   return nlhdlrdata->nexprs;
1556 	}
1557 	
1558 	/** adds a globally valid inequality of the form \f$\text{xcoef}\cdot x \leq \text{ycoef} \cdot y + \text{constant}\f$ to a product expression of the form \f$x\cdot y\f$ */
1559 	SCIP_RETCODE SCIPaddIneqBilinear(
1560 	   SCIP*                 scip,               /**< SCIP data structure */
1561 	   SCIP_NLHDLR*          nlhdlr,             /**< nonlinear handler */
1562 	   SCIP_EXPR*            expr,               /**< product expression */
1563 	   SCIP_Real             xcoef,              /**< x coefficient */
1564 	   SCIP_Real             ycoef,              /**< y coefficient */
1565 	   SCIP_Real             constant,           /**< constant part */
1566 	   SCIP_Bool*            success             /**< buffer to store whether inequality has been accepted */
1567 	   )
1568 	{
1569 	   SCIP_NLHDLREXPRDATA* nlhdlrexprdata;
1570 	   SCIP_VAR* x;
1571 	   SCIP_VAR* y;
1572 	   SCIP_Real* ineqs;
1573 	   SCIP_Real viol1;
1574 	   SCIP_Real viol2;
1575 	   SCIP_Bool underestimate;
1576 	   int nineqs;
1577 	   int i;
1578 	
1579 	   assert(scip != NULL);
1580 	   assert(nlhdlr != NULL);
1581 	   assert(strcmp(SCIPnlhdlrGetName(nlhdlr), NLHDLR_NAME) == 0);
1582 	   assert(expr != NULL);
1583 	   assert(SCIPexprGetNChildren(expr) == 2);
1584 	   assert(xcoef != SCIP_INVALID); /*lint !e777 */
1585 	   assert(ycoef != SCIP_INVALID); /*lint !e777 */
1586 	   assert(constant != SCIP_INVALID); /*lint !e777 */
1587 	   assert(success != NULL);
1588 	
1589 	   *success = FALSE;
1590 	
1591 	   /* find nonlinear handler expression handler data */
1592 	   nlhdlrexprdata = SCIPgetNlhdlrExprDataNonlinear(nlhdlr, expr);
1593 	
1594 	   if( nlhdlrexprdata == NULL )
1595 	   {
1596 	      SCIPwarningMessage(scip, "nonlinear expression data has not been found. "
1597 	            "Skip SCIPaddConsExprExprProductBilinearIneq()\n");
1598 	      return SCIP_OKAY;
1599 	   }
1600 	
1601 	   /* ignore inequalities that only yield to a (possible) bound tightening */
1602 	   if( SCIPisFeasZero(scip, xcoef) || SCIPisFeasZero(scip, ycoef) )
1603 	      return SCIP_OKAY;
1604 	
1605 	   /* collect variables */
1606 	   x = SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(expr)[0]);
1607 	   y = SCIPgetExprAuxVarNonlinear(SCIPexprGetChildren(expr)[1]);
1608 	   assert(x != NULL);
1609 	   assert(y != NULL);
1610 	   assert(x != y);
1611 	
1612 	   /* normalize inequality s.t. xcoef in {-1,1} */
1613 	   if( !SCIPisEQ(scip, REALABS(xcoef), 1.0) )
1614 	   {
1615 	      constant /= REALABS(xcoef);
1616 	      ycoef /= REALABS(xcoef);
1617 	      xcoef /= REALABS(xcoef);
1618 	   }
1619 	
1620 	   /* coefficients of the inequality determine whether the inequality can be used for under- or overestimation */
1621 	   underestimate = xcoef * ycoef > 0;
1622 	
1623 	   SCIPdebugMsg(scip, "add inequality for a bilinear term: %g %s <= %g %s + %g (underestimate=%u)\n", xcoef,
1624 	      SCIPvarGetName(x), ycoef, SCIPvarGetName(y), constant, underestimate);
1625 	
1626 	   /* compute violation of the inequality of the important corner points */
1627 	   getIneqViol(x, y, xcoef, ycoef, constant, &viol1, &viol2);
1628 	   SCIPdebugMsg(scip, "violations of inequality = (%g,%g)\n", viol1, viol2);
1629 	
1630 	   /* inequality does not cutoff one of the important corner points -> skip */
1631 	   if( SCIPisFeasLE(scip, MAX(viol1, viol2), 0.0) )
1632 	      return SCIP_OKAY;
1633 	
1634 	   if( underestimate )
1635 	   {
1636 	      ineqs = nlhdlrexprdata->underineqs;
1637 	      nineqs = nlhdlrexprdata->nunderineqs;
1638 	   }
1639 	   else
1640 	   {
1641 	      ineqs = nlhdlrexprdata->overineqs;
1642 	      nineqs = nlhdlrexprdata->noverineqs;
1643 	   }
1644 	   assert( nineqs >= 0 );
1645 	
1646 	   /* check for a duplicate */
1647 	   for( i = 0; i < nineqs; ++i )
1648 	   {
1649 	      if( SCIPisFeasEQ(scip, xcoef, ineqs[3*i]) && SCIPisFeasEQ(scip, ycoef, ineqs[3*i+1]) /*lint !e661*/
1650 	         && SCIPisFeasEQ(scip, constant, ineqs[3*i+2]) )
1651 	      {
1652 	         SCIPdebugMsg(scip, "inequality already found -> skip\n");
1653 	         return SCIP_OKAY;
1654 	      }
1655 	   }
1656 	
1657 	   {
1658 	      SCIP_Real viols1[2] = {0.0, 0.0};
1659 	      SCIP_Real viols2[2] = {0.0, 0.0};
1660 	
1661 	      /* compute violations of existing inequalities */
1662 	      for( i = 0; i < nineqs; ++i )
1663 	      {
1664 	         getIneqViol(x, y, ineqs[3*i], ineqs[3*i+1], ineqs[3*i+2], &viols1[i], &viols2[i]); /*lint !e661*/
1665 	
1666 	         /* check whether an existing inequality is dominating the candidate */
1667 	         if( SCIPisGE(scip, viols1[i], viol1) && SCIPisGE(scip, viols2[i], viol2) ) /*lint !e661*/
1668 	         {
1669 	            SCIPdebugMsg(scip, "inequality is dominated by %d -> skip\n", i);
1670 	            return SCIP_OKAY;
1671 	         }
1672 	
1673 	         /* replace inequality if candidate is dominating it */
1674 	         if( SCIPisLT(scip, viols1[i], viol1) && SCIPisLT(scip, viols2[i], viol2) ) /*lint !e661*/
1675 	         {
1676 	            SCIPdebugMsg(scip, "inequality dominates %d -> replace\n", i);
1677 	            ineqs[3*i] = xcoef; /*lint !e661*/
1678 	            ineqs[3*i+1] = ycoef; /*lint !e661*/
1679 	            ineqs[3*i+2] = constant; /*lint !e661*/
1680 	            *success = TRUE;
1681 	         }
1682 	      }
1683 	
1684 	      /* inequality is not dominated by other inequalities -> add if we have less than 2 inequalities */
1685 	      if( nineqs < 2 )
1686 	      {
1687 	         ineqs[3*nineqs] = xcoef;
1688 	         ineqs[3*nineqs + 1] = ycoef;
1689 	         ineqs[3*nineqs + 2] = constant;
1690 	         *success = TRUE;
1691 	         SCIPdebugMsg(scip, "add inequality\n");
1692 	
1693 	         /* increase number of inequalities */
1694 	         if( underestimate )
1695 	            ++(nlhdlrexprdata->nunderineqs);
1696 	         else
1697 	            ++(nlhdlrexprdata->noverineqs);
1698 	      }
1699 	   }
1700 	
1701 	   if( *success )
1702 	   {
1703 	      /* With the added inequalities, we can potentially compute tighter activities for the expression,
1704 	       * so constraints that contain this expression should be propagated again.
1705 	       * We don't have a direct expression to constraint mapping, though. This call marks all expr-constraints
1706 	       * which include any of the variables that this expression depends on for propagation.
1707 	       */
1708 	      SCIP_CALL( SCIPmarkExprPropagateNonlinear(scip, expr) );
1709 	   }
1710 	
1711 	   return SCIP_OKAY;
1712 	}
1713 	
1714 	/** computes coefficients of linearization of a bilinear term in a reference point */
1715 	void SCIPaddBilinLinearization(
1716 	   SCIP*                 scip,               /**< SCIP data structure */
1717 	   SCIP_Real             bilincoef,          /**< coefficient of bilinear term */
1718 	   SCIP_Real             refpointx,          /**< point where to linearize first  variable */
1719 	   SCIP_Real             refpointy,          /**< point where to linearize second variable */
1720 	   SCIP_Real*            lincoefx,           /**< buffer to add coefficient of first  variable in linearization */
1721 	   SCIP_Real*            lincoefy,           /**< buffer to add coefficient of second variable in linearization */
1722 	   SCIP_Real*            linconstant,        /**< buffer to add constant of linearization */
1723 	   SCIP_Bool*            success             /**< buffer to set to FALSE if linearization has failed due to large numbers */
1724 	   )
1725 	{
1726 	   SCIP_Real constant;
1727 	
1728 	   assert(scip != NULL);
1729 	   assert(lincoefx != NULL);
1730 	   assert(lincoefy != NULL);
1731 	   assert(linconstant != NULL);
1732 	   assert(success != NULL);
1733 	
1734 	   if( bilincoef == 0.0 )
1735 	      return;
1736 	
1737 	   if( SCIPisInfinity(scip, REALABS(refpointx)) || SCIPisInfinity(scip, REALABS(refpointy)) )
1738 	   {
1739 	      *success = FALSE;
1740 	      return;
1741 	   }
1742 	
1743 	   /* bilincoef * x * y ->  bilincoef * (refpointx * refpointy + refpointy * (x - refpointx) + refpointx * (y - refpointy))
1744 	    *                    = -bilincoef * refpointx * refpointy + bilincoef * refpointy * x + bilincoef * refpointx * y
1745 	    */
1746 	
1747 	   constant = -bilincoef * refpointx * refpointy;
1748 	
1749 	   if( SCIPisInfinity(scip, REALABS(bilincoef * refpointx)) || SCIPisInfinity(scip, REALABS(bilincoef * refpointy))
1750 	      || SCIPisInfinity(scip, REALABS(constant)) )
1751 	   {
1752 	      *success = FALSE;
1753 	      return;
1754 	   }
1755 	
1756 	   *lincoefx    += bilincoef * refpointy;
1757 	   *lincoefy    += bilincoef * refpointx;
1758 	   *linconstant += constant;
1759 	}
1760 	
1761 	/** computes coefficients of McCormick under- or overestimation of a bilinear term */
1762 	void SCIPaddBilinMcCormick(
1763 	   SCIP*                 scip,               /**< SCIP data structure */
1764 	   SCIP_Real             bilincoef,          /**< coefficient of bilinear term */
1765 	   SCIP_Real             lbx,                /**< lower bound on first variable */
1766 	   SCIP_Real             ubx,                /**< upper bound on first variable */
1767 	   SCIP_Real             refpointx,          /**< reference point for first variable */
1768 	   SCIP_Real             lby,                /**< lower bound on second variable */
1769 	   SCIP_Real             uby,                /**< upper bound on second variable */
1770 	   SCIP_Real             refpointy,          /**< reference point for second variable */
1771 	   SCIP_Bool             overestimate,       /**< whether to compute an overestimator instead of an underestimator */
1772 	   SCIP_Real*            lincoefx,           /**< buffer to add coefficient of first  variable in linearization */
1773 	   SCIP_Real*            lincoefy,           /**< buffer to add coefficient of second variable in linearization */
1774 	   SCIP_Real*            linconstant,        /**< buffer to add constant of linearization */
1775 	   SCIP_Bool*            success             /**< buffer to set to FALSE if linearization has failed due to large numbers */
1776 	   )
1777 	{
1778 	   SCIP_Real constant;
1779 	   SCIP_Real coefx;
1780 	   SCIP_Real coefy;
1781 	
1782 	   assert(scip != NULL);
1783 	   assert(!SCIPisInfinity(scip,  lbx));
1784 	   assert(!SCIPisInfinity(scip, -ubx));
1785 	   assert(!SCIPisInfinity(scip,  lby));
1786 	   assert(!SCIPisInfinity(scip, -uby));
1787 	   assert(SCIPisInfinity(scip,  -lbx) || SCIPisLE(scip, lbx, ubx));
1788 	   assert(SCIPisInfinity(scip,  -lby) || SCIPisLE(scip, lby, uby));
1789 	   assert(SCIPisInfinity(scip,  -lbx) || SCIPisLE(scip, lbx, refpointx));
1790 	   assert(SCIPisInfinity(scip,  -lby) || SCIPisLE(scip, lby, refpointy));
1791 	   assert(SCIPisInfinity(scip,  ubx) || SCIPisGE(scip, ubx, refpointx));
1792 	   assert(SCIPisInfinity(scip,  uby) || SCIPisGE(scip, uby, refpointy));
1793 	   assert(lincoefx != NULL);
1794 	   assert(lincoefy != NULL);
1795 	   assert(linconstant != NULL);
1796 	   assert(success != NULL);
1797 	
1798 	   if( bilincoef == 0.0 )
1799 	      return;
1800 	
1801 	   if( overestimate )
1802 	      bilincoef = -bilincoef;
1803 	
1804 	   if( SCIPisRelEQ(scip, lbx, ubx) && SCIPisRelEQ(scip, lby, uby) )
1805 	   {
1806 	      /* both x and y are mostly fixed */
1807 	      SCIP_Real cand1;
1808 	      SCIP_Real cand2;
1809 	      SCIP_Real cand3;
1810 	      SCIP_Real cand4;
1811 	
1812 	      coefx = 0.0;
1813 	      coefy = 0.0;
1814 	
1815 	      /* estimate x * y by constant */
1816 	      cand1 = lbx * lby;
1817 	      cand2 = lbx * uby;
1818 	      cand3 = ubx * lby;
1819 	      cand4 = ubx * uby;
1820 	
1821 	      /* take most conservative value for underestimator */
1822 	      if( bilincoef < 0.0 )
1823 	         constant = bilincoef * MAX( MAX(cand1, cand2), MAX(cand3, cand4) );
1824 	      else
1825 	         constant = bilincoef * MIN( MIN(cand1, cand2), MIN(cand3, cand4) );
1826 	   }
1827 	   else if( bilincoef > 0.0 )
1828 	   {
1829 	      /* either x or y is not fixed and coef > 0.0 */
1830 	      if( !SCIPisInfinity(scip, -lbx) && !SCIPisInfinity(scip, -lby) &&
1831 	         (SCIPisInfinity(scip,  ubx) || SCIPisInfinity(scip,  uby)
1832 	            || (uby - refpointy) * (ubx - refpointx) >= (refpointy - lby) * (refpointx - lbx)) )
1833 	      {
1834 	         if( SCIPisRelEQ(scip, lbx, ubx) )
1835 	         {
1836 	            /* x*y = lbx * y + (x-lbx) * y >= lbx * y + (x-lbx) * lby >= lbx * y + min{(ubx-lbx) * lby, 0 * lby} */
1837 	            coefx    =  0.0;
1838 	            coefy    =  bilincoef * lbx;
1839 	            constant =  bilincoef * (lby < 0.0 ? (ubx-lbx) * lby : 0.0);
1840 	         }
1841 	         else if( SCIPisRelEQ(scip, lby, uby) )
1842 	         {
1843 	            /* x*y = lby * x + (y-lby) * x >= lby * x + (y-lby) * lbx >= lby * x + min{(uby-lby) * lbx, 0 * lbx} */
1844 	            coefx    =  bilincoef * lby;
1845 	            coefy    =  0.0;
1846 	            constant =  bilincoef * (lbx < 0.0 ? (uby-lby) * lbx : 0.0);
1847 	         }
1848 	         else
1849 	         {
1850 	            coefx    =  bilincoef * lby;
1851 	            coefy    =  bilincoef * lbx;
1852 	            constant = -bilincoef * lbx * lby;
1853 	         }
1854 	      }
1855 	      else if( !SCIPisInfinity(scip, ubx) && !SCIPisInfinity(scip, uby) )
1856 	      {
1857 	         if( SCIPisRelEQ(scip, lbx, ubx) )
1858 	         {
1859 	            /* x*y = ubx * y + (x-ubx) * y >= ubx * y + (x-ubx) * uby >= ubx * y + min{(lbx-ubx) * uby, 0 * uby} */
1860 	            coefx    =  0.0;
1861 	            coefy    =  bilincoef * ubx;
1862 	            constant =  bilincoef * (uby > 0.0 ? (lbx - ubx) * uby : 0.0);
1863 	         }
1864 	         else if( SCIPisRelEQ(scip, lby, uby) )
1865 	         {
1866 	            /* x*y = uby * x + (y-uby) * x >= uby * x + (y-uby) * ubx >= uby * x + min{(lby-uby) * ubx, 0 * ubx} */
1867 	            coefx    =  bilincoef * uby;
1868 	            coefy    =  0.0;
1869 	            constant =  bilincoef * (ubx > 0.0 ? (lby - uby) * ubx : 0.0);
1870 	         }
1871 	         else
1872 	         {
1873 	            coefx    =  bilincoef * uby;
1874 	            coefy    =  bilincoef * ubx;
1875 	            constant = -bilincoef * ubx * uby;
1876 	         }
1877 	      }
1878 	      else
1879 	      {
1880 	         *success = FALSE;
1881 	         return;
1882 	      }
1883 	   }
1884 	   else
1885 	   {
1886 	      /* either x or y is not fixed and coef < 0.0 */
1887 	      if( !SCIPisInfinity(scip,  ubx) && !SCIPisInfinity(scip, -lby) &&
1888 	         (SCIPisInfinity(scip, -lbx) || SCIPisInfinity(scip,  uby)
1889 	            || (ubx - lbx) * (refpointy - lby) <= (uby - lby) * (refpointx - lbx)) )
1890 	      {
1891 	         if( SCIPisRelEQ(scip, lbx, ubx) )
1892 	         {
1893 	            /* x*y = ubx * y + (x-ubx) * y <= ubx * y + (x-ubx) * lby <= ubx * y + max{(lbx-ubx) * lby, 0 * lby} */
1894 	            coefx    =  0.0;
1895 	            coefy    =  bilincoef * ubx;
1896 	            constant =  bilincoef * (lby < 0.0 ? (lbx - ubx) * lby : 0.0);
1897 	         }
1898 	         else if( SCIPisRelEQ(scip, lby, uby) )
1899 	         {
1900 	            /* x*y = lby * x + (y-lby) * x <= lby * x + (y-lby) * ubx <= lby * x + max{(uby-lby) * ubx, 0 * ubx} */
1901 	            coefx    =  bilincoef * lby;
1902 	            coefy    =  0.0;
1903 	            constant =  bilincoef * (ubx > 0.0 ? (uby - lby) * ubx : 0.0);
1904 	         }
1905 	         else
1906 	         {
1907 	            coefx    =  bilincoef * lby;
1908 	            coefy    =  bilincoef * ubx;
1909 	            constant = -bilincoef * ubx * lby;
1910 	         }
1911 	      }
1912 	      else if( !SCIPisInfinity(scip, -lbx) && !SCIPisInfinity(scip, uby) )
1913 	      {
1914 	         if( SCIPisRelEQ(scip, lbx, ubx) )
1915 	         {
1916 	            /* x*y = lbx * y + (x-lbx) * y <= lbx * y + (x-lbx) * uby <= lbx * y + max{(ubx-lbx) * uby, 0 * uby} */
1917 	            coefx    =  0.0;
1918 	            coefy    =  bilincoef * lbx;
1919 	            constant =  bilincoef * (uby > 0.0 ? (ubx - lbx) * uby : 0.0);
1920 	         }
1921 	         else if( SCIPisRelEQ(scip, lby, uby) )
1922 	         {
1923 	            /* x*y = uby * x + (y-uby) * x <= uby * x + (y-uby) * lbx <= uby * x + max{(lby-uby) * lbx, 0 * lbx} */
1924 	            coefx    =  bilincoef * uby;
1925 	            coefy    =  0.0;
1926 	            constant =  bilincoef * (lbx < 0.0 ? (lby - uby) * lbx : 0.0);
1927 	         }
1928 	         else
1929 	         {
1930 	            coefx    =  bilincoef * uby;
1931 	            coefy    =  bilincoef * lbx;
1932 	            constant = -bilincoef * lbx * uby;
1933 	         }
1934 	      }
1935 	      else
1936 	      {
1937 	         *success = FALSE;
1938 	         return;
1939 	      }
1940 	   }
1941 	
1942 	   if( SCIPisInfinity(scip, REALABS(coefx)) || SCIPisInfinity(scip, REALABS(coefy))
1943 	      || SCIPisInfinity(scip, REALABS(constant)) )
1944 	   {
1945 	      *success = FALSE;
1946 	      return;
1947 	   }
1948 	
1949 	   if( overestimate )
1950 	   {
1951 	      coefx    = -coefx;
1952 	      coefy    = -coefy;
1953 	      constant = -constant;
1954 	   }
1955 	
1956 	   SCIPdebugMsg(scip, "%.15g * x[%.15g,%.15g] * y[%.15g,%.15g] %c= %.15g * x %+.15g * y %+.15g\n", bilincoef, lbx, ubx,
1957 	      lby, uby, overestimate ? '<' : '>', coefx, coefy, constant);
1958 	
1959 	   *lincoefx    += coefx;
1960 	   *lincoefy    += coefy;
1961 	   *linconstant += constant;
1962 	}
1963 	
1964 	/** computes coefficients of linearization of a bilinear term in a reference point when given a linear inequality
1965 	 *  involving only the variables of the bilinear term
1966 	 *
1967 	 *  @note the formulas are extracted from "Convex envelopes of bivariate functions through the solution of KKT systems"
1968 	 *        by Marco Locatelli
1969 	 */
1970 	void SCIPcomputeBilinEnvelope1(
1971 	   SCIP*                 scip,               /**< SCIP data structure */
1972 	   SCIP_Real             bilincoef,          /**< coefficient of bilinear term */
1973 	   SCIP_Real             lbx,                /**< lower bound on first variable */
1974 	   SCIP_Real             ubx,                /**< upper bound on first variable */
1975 	   SCIP_Real             refpointx,          /**< reference point for first variable */
1976 	   SCIP_Real             lby,                /**< lower bound on second variable */
1977 	   SCIP_Real             uby,                /**< upper bound on second variable */
1978 	   SCIP_Real             refpointy,          /**< reference point for second variable */
1979 	   SCIP_Bool             overestimate,       /**< whether to compute an overestimator instead of an underestimator */
1980 	   SCIP_Real             xcoef,              /**< x coefficient of linear inequality; must be in {-1,0,1} */
1981 	   SCIP_Real             ycoef,              /**< y coefficient of linear inequality */
1982 	   SCIP_Real             constant,           /**< constant of linear inequality */
1983 	   SCIP_Real* RESTRICT   lincoefx,           /**< buffer to store coefficient of first  variable in linearization */
1984 	   SCIP_Real* RESTRICT   lincoefy,           /**< buffer to store coefficient of second variable in linearization */
1985 	   SCIP_Real* RESTRICT   linconstant,        /**< buffer to store constant of linearization */
1986 	   SCIP_Bool* RESTRICT   success             /**< buffer to store whether linearization was successful */
1987 	   )
1988 	{
1989 	   SCIP_Real xs[2] = {lbx, ubx};
1990 	   SCIP_Real ys[2] = {lby, uby};
1991 	   SCIP_Real minx;
1992 	   SCIP_Real maxx;
1993 	   SCIP_Real miny;
1994 	   SCIP_Real maxy;
1995 	   SCIP_Real QUAD(lincoefyq);
1996 	   SCIP_Real QUAD(lincoefxq);
1997 	   SCIP_Real QUAD(linconstantq);
1998 	   SCIP_Real QUAD(denomq);
1999 	   SCIP_Real QUAD(mjq);
2000 	   SCIP_Real QUAD(qjq);
2001 	   SCIP_Real QUAD(xjq);
2002 	   SCIP_Real QUAD(yjq);
2003 	   SCIP_Real QUAD(tmpq);
2004 	   SCIP_Real vx;
2005 	   SCIP_Real vy;
2006 	   int n;
2007 	   int i;
2008 	
2009 	   assert(scip != NULL);
2010 	   assert(!SCIPisInfinity(scip,  lbx));
2011 	   assert(!SCIPisInfinity(scip, -ubx));
2012 	   assert(!SCIPisInfinity(scip,  lby));
2013 	   assert(!SCIPisInfinity(scip, -uby));
2014 	   assert(SCIPisLE(scip, lbx, ubx));
2015 	   assert(SCIPisLE(scip, lby, uby));
2016 	   assert(SCIPisLE(scip, lbx, refpointx));
2017 	   assert(SCIPisGE(scip, ubx, refpointx));
2018 	   assert(SCIPisLE(scip, lby, refpointy));
2019 	   assert(SCIPisGE(scip, uby, refpointy));
2020 	   assert(lincoefx != NULL);
2021 	   assert(lincoefy != NULL);
2022 	   assert(linconstant != NULL);
2023 	   assert(success != NULL);
2024 	   assert(xcoef == 0.0 || xcoef == -1.0 || xcoef == 1.0); /*lint !e777*/
2025 	   assert(ycoef != SCIP_INVALID && ycoef != 0.0); /*lint !e777*/
2026 	   assert(constant != SCIP_INVALID); /*lint !e777*/
2027 	
2028 	   *success = FALSE;
2029 	   *lincoefx = SCIP_INVALID;
2030 	   *lincoefy = SCIP_INVALID;
2031 	   *linconstant = SCIP_INVALID;
2032 	
2033 	   /* reference point does not satisfy linear inequality */
2034 	   if( SCIPisFeasGT(scip, xcoef * refpointx - ycoef * refpointy - constant, 0.0) )
2035 	      return;
2036 	
2037 	   /* compute minimal and maximal bounds on x and y for accepting the reference point */
2038 	   minx = lbx + 0.01 * (ubx-lbx);
2039 	   maxx = ubx - 0.01 * (ubx-lbx);
2040 	   miny = lby + 0.01 * (uby-lby);
2041 	   maxy = uby - 0.01 * (uby-lby);
2042 	
2043 	   /* check whether the reference point is in [minx,maxx]x[miny,maxy] */
2044 	   if( SCIPisLE(scip, refpointx, minx) || SCIPisGE(scip, refpointx, maxx)
2045 	      || SCIPisLE(scip, refpointy, miny) || SCIPisGE(scip, refpointy, maxy) )
2046 	      return;
2047 	
2048 	   /* always consider xy without the bilinear coefficient */
2049 	   if( bilincoef < 0.0 )
2050 	      overestimate = !overestimate;
2051 	
2052 	   /* we use same notation as in "Convex envelopes of bivariate functions through the solution of KKT systems", 2016 */
2053 	   /* mj = xcoef / ycoef */
2054 	   SCIPquadprecDivDD(mjq, xcoef, ycoef);
2055 	
2056 	   /* qj = -constant / ycoef */
2057 	   SCIPquadprecDivDD(qjq, -constant, ycoef);
2058 	
2059 	   /* mj > 0 => underestimate; mj < 0 => overestimate */
2060 	   if( SCIPisNegative(scip, QUAD_TO_DBL(mjq)) != overestimate )
2061 	      return;
2062 	
2063 	   /* get the corner point that satisfies the linear inequality xcoef*x <= ycoef*y + constant */
2064 	   if( !overestimate )
2065 	   {
2066 	      ys[0] = uby;
2067 	      ys[1] = lby;
2068 	   }
2069 	
2070 	   vx = SCIP_INVALID;
2071 	   vy = SCIP_INVALID;
2072 	   n = 0;
2073 	   for( i = 0; i < 2; ++i )
2074 	   {
2075 	      SCIP_Real activity = xcoef * xs[i] - ycoef * ys[i] - constant;
2076 	      if( SCIPisLE(scip, activity, 0.0) )
2077 	      {
2078 	         /* corner point is satisfies inequality */
2079 	         vx = xs[i];
2080 	         vy = ys[i];
2081 	      }
2082 	      else if( SCIPisFeasGT(scip, activity, 0.0) )
2083 	         /* corner point is clearly cut off */
2084 	         ++n;
2085 	   }
2086 	
2087 	   /* skip if no corner point satisfies the inequality or if no corner point is cut off
2088 	    * (that is, all corner points satisfy the inequality almost [1e-9..1e-6]) */
2089 	   if( n != 1 || vx == SCIP_INVALID || vy == SCIP_INVALID ) /*lint !e777*/
2090 	      return;
2091 	
2092 	   /* denom = mj*(refpointx - vx) + vy - refpointy */
2093 	   SCIPquadprecSumDD(denomq, refpointx, -vx); /* refpoint - vx */
2094 	   SCIPquadprecProdQQ(denomq, denomq, mjq); /* mj * (refpoint - vx) */
2095 	   SCIPquadprecSumQD(denomq, denomq, vy); /* mj * (refpoint - vx) + vy */
2096 	   SCIPquadprecSumQD(denomq, denomq, -refpointy); /* mj * (refpoint - vx) + vy - refpointy */
2097 	
2098 	   if( SCIPisZero(scip, QUAD_TO_DBL(denomq)) )
2099 	      return;
2100 	
2101 	   /* (xj,yj) is the projection onto the line xcoef*x = ycoef*y + constant */
2102 	   /* xj = (refpointx*(vy - qj) - vx*(refpointy - qj)) / denom */
2103 	   SCIPquadprecProdQD(xjq, qjq, -1.0); /* - qj */
2104 	   SCIPquadprecSumQD(xjq, xjq, vy); /* vy - qj */
2105 	   SCIPquadprecProdQD(xjq, xjq, refpointx); /* refpointx * (vy - qj) */
2106 	   SCIPquadprecProdQD(tmpq, qjq, -1.0); /* - qj */
2107 	   SCIPquadprecSumQD(tmpq, tmpq, refpointy); /* refpointy - qj */
2108 	   SCIPquadprecProdQD(tmpq, tmpq, -vx); /* - vx * (refpointy - qj) */
2109 	   SCIPquadprecSumQQ(xjq, xjq, tmpq); /* refpointx * (vy - qj) - vx * (refpointy - qj) */
2110 	   SCIPquadprecDivQQ(xjq, xjq, denomq); /* (refpointx * (vy - qj) - vx * (refpointy - qj)) / denom */
2111 	
2112 	   /* yj = mj * xj + qj */
2113 	   SCIPquadprecProdQQ(yjq, mjq, xjq);
2114 	   SCIPquadprecSumQQ(yjq, yjq, qjq);
2115 	
2116 	   assert(SCIPisFeasEQ(scip, xcoef*QUAD_TO_DBL(xjq) - ycoef*QUAD_TO_DBL(yjq) - constant, 0.0));
2117 	
2118 	   /* check whether the projection is in [minx,maxx] x [miny,maxy]; this avoids numerical difficulties when the
2119 	    * projection is close to the variable bounds
2120 	    */
2121 	   if( SCIPisLE(scip, QUAD_TO_DBL(xjq), minx) || SCIPisGE(scip, QUAD_TO_DBL(xjq), maxx)
2122 	      || SCIPisLE(scip, QUAD_TO_DBL(yjq), miny) || SCIPisGE(scip, QUAD_TO_DBL(yjq), maxy) )
2123 	      return;
2124 	
2125 	   assert(vy - QUAD_TO_DBL(mjq)*vx - QUAD_TO_DBL(qjq) != 0.0);
2126 	
2127 	   /* lincoefy = (mj*SQR(xj) - 2.0*mj*vx*xj - qj*vx + vx*vy) / (vy - mj*vx - qj) */
2128 	   SCIPquadprecSquareQ(lincoefyq, xjq); /* xj^2 */
2129 	   SCIPquadprecProdQQ(lincoefyq, lincoefyq, mjq); /* mj * xj^2 */
2130 	   SCIPquadprecProdQQ(tmpq, mjq, xjq); /* mj * xj */
2131 	   SCIPquadprecProdQD(tmpq, tmpq, -2.0 * vx); /* -2 * vx * mj * xj */
2132 	   SCIPquadprecSumQQ(lincoefyq, lincoefyq, tmpq); /* mj * xj^2 -2 * vx * mj * xj */
2133 	   SCIPquadprecProdQD(tmpq, qjq, -vx); /* -qj * vx */
2134 	   SCIPquadprecSumQQ(lincoefyq, lincoefyq, tmpq); /* mj * xj^2 -2 * vx * mj * xj -qj * vx */
2135 	   SCIPquadprecProdDD(tmpq, vx, vy); /* vx * vy */
2136 	   SCIPquadprecSumQQ(lincoefyq, lincoefyq, tmpq); /* mj * xj^2 -2 * vx * mj * xj -qj * vx + vx * vy */
2137 	   SCIPquadprecProdQD(tmpq, mjq, vx); /* mj * vx */
2138 	   SCIPquadprecSumQD(tmpq, tmpq, -vy); /* -vy + mj * vx */
2139 	   SCIPquadprecSumQQ(tmpq, tmpq, qjq); /* -vy + mj * vx + qj */
2140 	   QUAD_SCALE(tmpq, -1.0); /* vy - mj * vx - qj */
2141 	   SCIPquadprecDivQQ(lincoefyq, lincoefyq, tmpq); /* (mj * xj^2 -2 * vx * mj * xj -qj * vx + vx * vy) / (vy - mj * vx - qj) */
2142 	
2143 	   /* lincoefx = 2.0*mj*xj + qj - mj*(*lincoefy) */
2144 	   SCIPquadprecProdQQ(lincoefxq, mjq, xjq); /* mj * xj */
2145 	   QUAD_SCALE(lincoefxq, 2.0); /* 2 * mj * xj */
2146 	   SCIPquadprecSumQQ(lincoefxq, lincoefxq, qjq); /* 2 * mj * xj + qj */
2147 	   SCIPquadprecProdQQ(tmpq, mjq, lincoefyq); /* mj * lincoefy */
2148 	   QUAD_SCALE(tmpq, -1.0); /* - mj * lincoefy */
2149 	   SCIPquadprecSumQQ(lincoefxq, lincoefxq, tmpq); /* 2 * mj * xj + qj - mj * lincoefy */
2150 	
2151 	   /* linconstant = -mj*SQR(xj) - (*lincoefy)*qj */
2152 	   SCIPquadprecSquareQ(linconstantq, xjq); /* xj^2 */
2153 	   SCIPquadprecProdQQ(linconstantq, linconstantq, mjq); /* mj * xj^2 */
2154 	   QUAD_SCALE(linconstantq, -1.0); /* - mj * xj^2 */
2155 	   SCIPquadprecProdQQ(tmpq, lincoefyq, qjq); /* lincoefy * qj */
2156 	   QUAD_SCALE(tmpq, -1.0); /* - lincoefy * qj */
2157 	   SCIPquadprecSumQQ(linconstantq, linconstantq, tmpq); /* - mj * xj^2 - lincoefy * qj */
2158 	
2159 	   /* consider the bilinear coefficient */
2160 	   SCIPquadprecProdQD(lincoefxq, lincoefxq, bilincoef);
2161 	   SCIPquadprecProdQD(lincoefyq, lincoefyq, bilincoef);
2162 	   SCIPquadprecProdQD(linconstantq, linconstantq, bilincoef);
2163 	   *lincoefx = QUAD_TO_DBL(lincoefxq);
2164 	   *lincoefy = QUAD_TO_DBL(lincoefyq);
2165 	   *linconstant = QUAD_TO_DBL(linconstantq);
2166 	
2167 	   /* cut needs to be tight at (vx,vy) and (xj,yj); otherwise we consider the cut to be numerically bad */
2168 	   *success = SCIPisFeasEQ(scip, (*lincoefx)*vx + (*lincoefy)*vy + (*linconstant), bilincoef*vx*vy)
2169 	      && SCIPisFeasEQ(scip, (*lincoefx)*QUAD_TO_DBL(xjq) + (*lincoefy)*QUAD_TO_DBL(yjq) + (*linconstant),
2170 	      bilincoef*QUAD_TO_DBL(xjq)*QUAD_TO_DBL(yjq));
2171 	
2172 	#ifndef NDEBUG
2173 	   {
2174 	      SCIP_Real activity = (*lincoefx)*refpointx + (*lincoefy)*refpointy + (*linconstant);
2175 	
2176 	      /* cut needs to under- or overestimate the bilinear term at the reference point */
2177 	      if( bilincoef < 0.0 )
2178 	         overestimate = !overestimate;
2179 	
2180 	      if( overestimate )
2181 	         assert(SCIPisFeasGE(scip, activity, bilincoef*refpointx*refpointy));
2182 	      else
2183 	         assert(SCIPisFeasLE(scip, activity, bilincoef*refpointx*refpointy));
2184 	   }
2185 	#endif
2186 	}
2187 	
2188 	/** computes coefficients of linearization of a bilinear term in a reference point when given two linear inequalities
2189 	 *  involving only the variables of the bilinear term
2190 	 *
2191 	 *  @note the formulas are extracted from "Convex envelopes of bivariate functions through the solution of KKT systems"
2192 	 *        by Marco Locatelli
2193 	 */
2194 	void SCIPcomputeBilinEnvelope2(
2195 	   SCIP*                 scip,               /**< SCIP data structure */
2196 	   SCIP_Real             bilincoef,          /**< coefficient of bilinear term */
2197 	   SCIP_Real             lbx,                /**< lower bound on first variable */
2198 	   SCIP_Real             ubx,                /**< upper bound on first variable */
2199 	   SCIP_Real             refpointx,          /**< reference point for first variable */
2200 	   SCIP_Real             lby,                /**< lower bound on second variable */
2201 	   SCIP_Real             uby,                /**< upper bound on second variable */
2202 	   SCIP_Real             refpointy,          /**< reference point for second variable */
2203 	   SCIP_Bool             overestimate,       /**< whether to compute an overestimator instead of an underestimator */
2204 	   SCIP_Real             xcoef1,             /**< x coefficient of linear inequality; must be in {-1,0,1} */
2205 	   SCIP_Real             ycoef1,             /**< y coefficient of linear inequality */
2206 	   SCIP_Real             constant1,          /**< constant of linear inequality */
2207 	   SCIP_Real             xcoef2,             /**< x coefficient of linear inequality; must be in {-1,0,1} */
2208 	   SCIP_Real             ycoef2,             /**< y coefficient of linear inequality */
2209 	   SCIP_Real             constant2,          /**< constant of linear inequality */
2210 	   SCIP_Real* RESTRICT   lincoefx,           /**< buffer to store coefficient of first  variable in linearization */
2211 	   SCIP_Real* RESTRICT   lincoefy,           /**< buffer to store coefficient of second variable in linearization */
2212 	   SCIP_Real* RESTRICT   linconstant,        /**< buffer to store constant of linearization */
2213 	   SCIP_Bool* RESTRICT   success             /**< buffer to store whether linearization was successful */
2214 	   )
2215 	{
2216 	   SCIP_Real mi, mj, qi, qj, xi, xj, yi, yj;
2217 	   SCIP_Real xcoef, ycoef, constant;
2218 	   SCIP_Real minx, maxx, miny, maxy;
2219 	
2220 	   assert(scip != NULL);
2221 	   assert(!SCIPisInfinity(scip,  lbx));
2222 	   assert(!SCIPisInfinity(scip, -ubx));
2223 	   assert(!SCIPisInfinity(scip,  lby));
2224 	   assert(!SCIPisInfinity(scip, -uby));
2225 	   assert(SCIPisLE(scip, lbx, ubx));
2226 	   assert(SCIPisLE(scip, lby, uby));
2227 	   assert(SCIPisLE(scip, lbx, refpointx));
2228 	   assert(SCIPisGE(scip, ubx, refpointx));
2229 	   assert(SCIPisLE(scip, lby, refpointy));
2230 	   assert(SCIPisGE(scip, uby, refpointy));
2231 	   assert(lincoefx != NULL);
2232 	   assert(lincoefy != NULL);
2233 	   assert(linconstant != NULL);
2234 	   assert(success != NULL);
2235 	   assert(xcoef1 != 0.0 && xcoef1 != SCIP_INVALID); /*lint !e777*/
2236 	   assert(ycoef1 != SCIP_INVALID && ycoef1 != 0.0); /*lint !e777*/
2237 	   assert(constant1 != SCIP_INVALID); /*lint !e777*/
2238 	   assert(xcoef2 != 0.0 && xcoef2 != SCIP_INVALID); /*lint !e777*/
2239 	   assert(ycoef2 != SCIP_INVALID && ycoef2 != 0.0); /*lint !e777*/
2240 	   assert(constant2 != SCIP_INVALID); /*lint !e777*/
2241 	
2242 	   *success = FALSE;
2243 	   *lincoefx = SCIP_INVALID;
2244 	   *lincoefy = SCIP_INVALID;
2245 	   *linconstant = SCIP_INVALID;
2246 	
2247 	   /* reference point does not satisfy linear inequalities */
2248 	   if( SCIPisFeasGT(scip, xcoef1 * refpointx - ycoef1 * refpointy - constant1, 0.0)
2249 	      || SCIPisFeasGT(scip, xcoef2 * refpointx - ycoef2 * refpointy - constant2, 0.0) )
2250 	      return;
2251 	
2252 	   /* compute minimal and maximal bounds on x and y for accepting the reference point */
2253 	   minx = lbx + 0.01 * (ubx-lbx);
2254 	   maxx = ubx - 0.01 * (ubx-lbx);
2255 	   miny = lby + 0.01 * (uby-lby);
2256 	   maxy = uby - 0.01 * (uby-lby);
2257 	
2258 	   /* check the reference point is in the interior of the domain */
2259 	   if( SCIPisLE(scip, refpointx, minx) || SCIPisGE(scip, refpointx, maxx)
2260 	      || SCIPisLE(scip, refpointy, miny) || SCIPisFeasGE(scip, refpointy, maxy) )
2261 	      return;
2262 	
2263 	   /* the sign of the x-coefficients of the two inequalities must be different; otherwise the convex or concave
2264 	    * envelope can be computed via SCIPcomputeBilinEnvelope1 for each inequality separately
2265 	    */
2266 	   if( (xcoef1 > 0) == (xcoef2 > 0) )
2267 	      return;
2268 	
2269 	   /* always consider xy without the bilinear coefficient */
2270 	   if( bilincoef < 0.0 )
2271 	      overestimate = !overestimate;
2272 	
2273 	   /* we use same notation as in "Convex envelopes of bivariate functions through the solution of KKT systems", 2016 */
2274 	   mi = xcoef1 / ycoef1;
2275 	   qi = -constant1 / ycoef1;
2276 	   mj = xcoef2 / ycoef2;
2277 	   qj = -constant2 / ycoef2;
2278 	
2279 	   /* mi, mj > 0 => underestimate; mi, mj < 0 => overestimate */
2280 	   if( SCIPisNegative(scip, mi) != overestimate || SCIPisNegative(scip, mj) != overestimate )
2281 	      return;
2282 	
2283 	   /* compute cut according to Locatelli 2016 */
2284 	   computeBilinEnvelope2(scip, refpointx, refpointy, mi, qi, mj, qj, &xi, &yi, &xj, &yj, &xcoef, &ycoef, &constant);
2285 	   assert(SCIPisRelEQ(scip, mi*xi + qi, yi));
2286 	   assert(SCIPisRelEQ(scip, mj*xj + qj, yj));
2287 	
2288 	   /* it might happen that (xi,yi) = (xj,yj) if the two lines intersect */
2289 	   if( SCIPisEQ(scip, xi, xj) && SCIPisEQ(scip, yi, yj) )
2290 	      return;
2291 	
2292 	   /* check whether projected points are in the interior */
2293 	   if( SCIPisLE(scip, xi, minx) || SCIPisGE(scip, xi, maxx) || SCIPisLE(scip, yi, miny) || SCIPisGE(scip, yi, maxy) )
2294 	      return;
2295 	   if( SCIPisLE(scip, xj, minx) || SCIPisGE(scip, xj, maxx) || SCIPisLE(scip, yj, miny) || SCIPisGE(scip, yj, maxy) )
2296 	      return;
2297 	
2298 	   *lincoefx = bilincoef * xcoef;
2299 	   *lincoefy = bilincoef * ycoef;
2300 	   *linconstant = bilincoef * constant;
2301 	
2302 	   /* cut needs to be tight at (vx,vy) and (xj,yj) */
2303 	   *success = SCIPisFeasEQ(scip, (*lincoefx)*xi + (*lincoefy)*yi + (*linconstant), bilincoef*xi*yi)
2304 	      && SCIPisFeasEQ(scip, (*lincoefx)*xj + (*lincoefy)*yj + (*linconstant), bilincoef*xj*yj);
2305 	
2306 	#ifndef NDEBUG
2307 	   {
2308 	      SCIP_Real activity = (*lincoefx)*refpointx + (*lincoefy)*refpointy + (*linconstant);
2309 	
2310 	      /* cut needs to under- or overestimate the bilinear term at the reference point */
2311 	      if( bilincoef < 0.0 )
2312 	         overestimate = !overestimate;
2313 	
2314 	      if( overestimate )
2315 	         assert(SCIPisFeasGE(scip, activity, bilincoef*refpointx*refpointy));
2316 	      else
2317 	         assert(SCIPisFeasLE(scip, activity, bilincoef*refpointx*refpointy));
2318 	   }
2319 	#endif
2320 	}
2321