3dbpp.cpp

/* ======================================================================
 3D BIN PACKING, Silvano Martello, David Pisinger, Daniele Vigo
 1998, 2003, 2006
 ====================================================================== */

/* This code solves the three-dimensional bin-packing problem, which
 * asks for an orthogonal packing of a given set of rectangular-shaped
 * boxes into the minimum number of three-dimensional rectangular bins.
 * Each box j=1,..,n is characterized by a width w_j, height h_j, and
 * depth d_j. An unlimited number of indentical three-dimensional bins,
 * having width W, height H and depth D is available. The boxes have fixed
 * orientation, i.e., they may not be rotated.
 *
 * A description of this code is found in the following papers:
 *
 *   S. Martello, D. Pisinger, D. Vigo, E. den Boef, J. Korst (2003)
 *   "Algorithms for General and Robot-packable Variants of the 
 *    Three-Dimensional Bin Packing Problem"
 *   submitted TOMS.
 * 
 *   S.Martello, D.Pisinger, D.Vigo (2000)
 *   "The three-dimensional bin packing problem"
 *   Operations Research, 48, 256-267
 *
 * The present code is written in ANSI-C, and has been compiled with
 * the GNU-C compiler using option "-ansi -pedantic" as well as the
 * HP-UX C compiler using option "-Aa" (ansi standard).
 *
 * This file contains the callable routine binpack3d with prototype
 *
 *   void binpack3d(int n, int W, int H, int D,
 *                  int *w, int *h, int *d, 
 *                  int *x, int *y, int *z, int *bno,
 *                  int *lb, int *ub, 
 *                  int nodelimit, int iterlimit, int timelimit, 
 *                  int *nodeused, int *iterused, int *timeused, 
 *                  int packingtype);
 *
 * the meaning of the parameters is the following:
 *   n         Size of problem, i.e., number of boxes to be packed.
 *             This value must be smaller than MAXBOXES defined below.
 *   W,H,D     Width, height and depth of every bin.
 *   w,h,d     Integer arrays of length n, where w[j], h[j], d[j]
 *             are the dimensions of box j for j=0,..,n-1.
 *   x,y,z,bno Integer arrays of length n where the solution found
 *             is returned. For each box j=0,..,n-1, the bin number
 *             it is packed into is given by bno[j], and x[j], y[j], z[j] 
 *             are the coordinates of it lower-left-backward corner.
 *   lb        Lower bound on the solution value (returned by the procedure).
 *   ub        Objective value of the solution found, i.e., number of bins
 *             used to pack the n boxes. (returned by the procedure).
 *   nodelimit maximum number of decision nodes to be explored in the
 *             main branching tree. If set to zero, the algorithm will
 *             run until an optimal solution is found (or timelimit or
 *             iterlimit is reached). Measured in thousands (see IUNIT).
 *   iterlimit maximum number of iterations in the ONEBIN algorithm
 *             which packs a single bin. If set to zero, the algorithm will
 *             run until an optimal solution is found (or timelimit or
 *             nodelimit is reached). Measured in thousands (see IUNIT).
 *   timelimit Time limit for solving the problem expressed in seconds.
 *             If set to zero, the algorithm will run until an optimal
 *             solution is found; otherwise it terminates after timelimit
 *             seconds with a heuristic solution. 
 *   nodeused  returns the number of branch-and-bound nodes investigated,
 *             measured in thousands (see IUNIT).
 *   iterused  returns the number of iterations in ONEBIN algorithm,
 *             measured in thousands (see IUNIT).
 *   timeused  returns the time used in miliseconds
 *   packingtype 
 *             Desired packing type. If set to zero, the algorithm will
 *             search for an optimal general packing; if set to one, it
 *             will search for a robot packing.
 *
 * (c) Copyright 1998, 2003, 2005
 *
 *   David Pisinger                        Silvano Martello, Daniele Vigo
 *   DIKU, University of Copenhagen        DEIS, University of Bologna
 *   Universitetsparken 1                  Viale Risorgimento 2
 *   Copenhagen, Denmark                   Bologna, Italy
 *
 * This code can be used free of charge for research and academic purposes 
 * only. 
 */

#include <stdlib.h>
#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <math.h>
#include <malloc.h>
#include <time.h>  
#include <limits.h>

#include <map>
#include "ext/xml_parser.h"
#include "ext/packlist.h"
#include "ext/response.h"

#include "3dbpp.h"

/* ======================================================================
 macros
 ====================================================================== */

#define TRUE  1           /* logical variables */
#define FALSE 0

#define WDIM  0           /* rotations of boxes */
#define HDIM  1
#define DDIM  2

#define GENERAL     0     /* packing type */
#define ROBOT       1

#define LEFT   0          /* relative placements */
#define RIGHT  1
#define UNDER  2
#define ABOVE  3
#define FRONT  4
#define BEHIND 5
#define UNDEF  6
#define RELMAX 8

#define STACKDEPTH (MAXBOXES*MAXBOXES*RELMAX)

#define VOL(i)            ((i)->w * (ptype) (i)->h * (i)->d)
#define MINIMUM(i,j)      ((i) < (j) ? (i) : (j))
#define MAXIMUM(i,j)      ((i) > (j) ? (i) : (j))
#define DIF(i,j)          ((int) ((j) - (i) + 1))
#define SWAPINT(a,b)      { register int t; t=*(a);*(a)=*(b);*(b)=t; }
#define SWAP(a,b)         { register box t; t=*(a);*(a)=*(b);*(b)=t; }
#define SWAPI(a,b)        { register itype t; t=(a);(a)=(b);(b)=t; }
#define SWAPP(a,b)        { register point t; t=*(a);*(a)=*(b);*(b)=t; }
#define DF(a,b)           ((r=(a).y-(b).y) != 0 ? r : (a).x-(b).x)

/* ======================================================================
 type declarations
 ====================================================================== */

typedef short boolean; /* logical variable      */
typedef short ntype; /* number of states,bins */
typedef short itype; /* can hold up to W,H,D  */
typedef long stype; /* can hold up to W*H*D  */
typedef long ptype; /* product multiplication */

/* box record */
typedef struct irec {
	ntype no; /* box number                            */
	itype w; /* box width  (x-size)                   */
	itype h; /* box height (y-size)                   */
	itype d; /* box depth  (z-size)                   */
	itype x; /* optimal x-position                    */
	itype y; /* optimal y-position                    */
	itype z; /* optimal z-position                    */
	ntype bno; /* bin number                            */
	boolean k; /* is the box chosen?                    */
	stype vol; /* volume of box                         */
	struct irec *ref; /* reference to original box (if necessary) */
} box;

/* all problem information */
typedef struct {
	itype W; /* x-size of bin                         */
	itype H; /* y-size of bin                         */
	itype D; /* z-size of bin                         */
	stype BVOL; /* volume of a bin                       */
	ntype n; /* number of boxes                       */
	boolean packtype; /* packing type: GENERAL or ROBOT        */
	box *fbox; /* first box in problem                  */
	box *lbox; /* last box in problem                   */
	box *fsol; /* first box in current solution         */
	box *lsol; /* last box in current solution          */
	box *fopt; /* first box in optimal solution         */
	box *lopt; /* last box in optimal solution          */
	boolean *closed; /* for each bin indicator whether closed */
	box *fclosed; /* first box in closed bins              */
	box *lclosed; /* last box in closed bins               */
	ntype noc; /* number of closed bins                 */
	itype mindim; /* currently smallest box length         */
	itype maxdim; /* currently largest box length          */
	stype maxfill; /* the best filling found                */
	int mcut; /* how many siblings at each node in b&b */

	/* different bounds */
	ntype bound0; /* Bound L_0 at root node                */
	ntype bound1; /* Bound L_1 at root node                */
	ntype bound2; /* Bound L_2 at root node                */
	ntype lb; /* best of the above                     */
	ntype z; /* currently best solution               */

	/* controle of 3d filler */
	int maxiter; /* max iterations in onebin_robot        */
	int miss; /* number boxes not packed in onebin_robot */

	/* debugging and controle information */
	int nodes; /* nodes in branch-and-bound             */
	int iterat; /* iterations in onebin_decision         */
	int subnodes; /* nodes in branch-and-bound             */
	int subiterat; /* iterations in onebin_decision         */
	int exfill; /* number of calls to onebin_decision    */
	int iter3d; /* iterations in onebin_robot or general */
	int zlayer; /* heuristic solution layer              */
	int zmcut; /* heuristic solution mcut               */
	double exacttopo; /* number of topological sorts           */
	double exacttopn; /* number of topological sorts           */
	int exactcall; /* number of calls to exact              */
	int exactn; /* largest problem for exact             */
	double genertime; /* time used in onebin_general           */
	double robottime; /* time used in onebin_robot             */
	double time; /* computing time                        */
	double lhtime; /* layer heuristic computing time        */
	double mhtime; /* mcut heuristic computing time         */
	int didpush; /* did the lower bound push up bound     */
	int maxclose; /* max number of closed bins at any time */
	int nodelimit; /* maximum number of nodes in main tree  */
	int iterlimit; /* maximum number of iterations in ONEBIN*/
	int timelimit; /* maximum amount of time to be used     */
} allinfo;

/* structure for greedy algorithm */
typedef struct {
	int lno; /* layer number                          */
	int d; /* depth of layer                        */
	int bno; /* bin no assigned to layer              */
	int z; /* z level of layer within bin           */
	int b; /* temporary bin number                  */
} heurpair;

/* structure for extreme points in a single bin */
typedef struct {
	itype x; /* x-coordinate                          */
	itype y; /* y-coordinate                          */
	itype z; /* z-coordinate                          */
} point;

/* structure for a domain pair in constraint programming */
typedef struct {
	int i; /* index of box i                        */
	int j; /* index of box j                        */
	int relation; /* relation between the two boxes        */
	boolean domain; /* domain of the two boxes               */
} domainpair;

/* set of domains */
typedef char domset[RELMAX];
typedef domset domline[MAXBOXES];

/* pointer to comparison function */
typedef int (*funcptr)(const void *, const void *);

/* ======================================================================
 global variables
 ====================================================================== */

/* boolean variable to indicate time-out situation */
boolean stopped;

/* counter used to ensure that 1D BPP at most performs MAXBPP iterations */
int bpiterat;

/* boolean variables to indicate when 1D packing algorithm should terminate */
boolean feasible, terminate;

/* stack of domain pairs */
domainpair domstack[STACKDEPTH];
domainpair *dompos, *domend;

/* domain of each box */
domline domain[MAXBOXES];

/* current relation between two boxes */
char relation[MAXBOXES][MAXBOXES];

/* debug variable to see level in recursive packing algorithm */
int bblevel;

/* =======================================================================
 error
 ======================================================================= */

void error(const char *str, ...) {
	va_list args;

	va_start(args, str);
	vprintf(str, args);
	printf("\n");
	va_end(args);
	printf("IRREGULAR PROGRAM TERMINATION\n");
	exit(-1);
}

/* **********************************************************************
 **********************************************************************
 Timing routines
 **********************************************************************
 ********************************************************************** */

/* This timing routine is based on the ANSI-C procedure "clock", which
 * has a resolution of 1000000 ticks per second. This however implies
 * that we pass the limit of a long integer after only 4295 seconds.
 * The following routine attempts to correct such situations by adding
 * the constant ULONG_MAX to the counter whenever wraparound can be
 * detected. But the user is advised to use a timing routine like "times"
 * (which however is not ANSI-C standard) for measuring longer time
 * periods.
 */

void timer(double *time) {
	static double tstart, tend, tprev;

	if (time == NULL) {
		clock(); /* one extra call to initialize clock */
		tstart = tprev = clock();
	} else {
		tend = clock();
		if (tend < tprev)
			tstart -= ULONG_MAX; /* wraparound occured */
		tprev = tend;
		*time = (tend - tstart) / CLOCKS_PER_SEC; /* convert to seconds */
	}
}

/* test for time limit */
void check_timelimit(long max) {
	double t;
	if (max == 0)
		return;
	timer(&t);
	if (t >= max) {
		if (!stopped)
			printf("TIMELIMIT\n");
		stopped = TRUE;
	}
}

/* test for node limit */
void check_nodelimit(long nodes, long max) {
	if (max == 0)
		return;
	if (nodes >= max) {
		if (!stopped)
			printf("NODELIMIT\n");
		stopped = TRUE;
	}
}

/* test for iteration limit */
void check_iterlimit(long iterations, long max) {
	if (max == 0)
		return;
	if (iterations >= max) {
		if (!stopped)
			printf("ITERLIMIT\n");
		stopped = TRUE;
	}
}

/* **********************************************************************
 **********************************************************************
 Small procedures
 **********************************************************************
 ********************************************************************** */

/* ======================================================================
 simple comparisions
 ====================================================================== */

/* Comparisons used as argument to qsort. */

int dcomp(box *a, box *b) {
	int r;
	r = b->d - a->d;
	if (r != 0)
		return r;
	else
		return b->no - a->no;
}
int hcomp(box *a, box *b) {
	int r;
	r = b->h - a->h;
	if (r != 0)
		return r;
	else
		return b->no - a->no;
}
int vcomp(box *a, box *b) /* volume decr. */
{
	int r;
	r = b->vol - a->vol;
	if (r != 0)
		return r;
	else
		return b->no - a->no;
}
int xcomp(heurpair *a, heurpair *b) /* depth decr. */
{
	int r;
	r = b->d - a->d;
	if (r != 0)
		return r;
	else
		return b->lno - a->lno;
}
int lcomp(heurpair *a, heurpair *b) /* layer number decr. */
{
	int r;
	r = a->lno - b->lno;
	if (r != 0)
		return r;
	else
		return b->d - a->d;
}

/* ======================================================================
 palloc
 ====================================================================== */

/* Memory allocation and freeing, with implicit check */

void *palloc(long sz, long no) {
	long size;
	void *p;

	size = sz * no;
	if (size == 0)
		size = 1;
	p = (void *) malloc(size);
	if (p == NULL)
		error("no memory size %ld", size);
	return p;
}

void pfree(void *p) {
	if (p == NULL)
		error("freeing null");
	free(p);
}

/* ======================================================================
 checksol
 ====================================================================== */

/* Check correctnes of solution, i.e., no boxes overlap, no duplicated boxes.
 * If the solution should be robot packable, it is checked that there
 * exists an ordering of the boxes such that they can be removed one
 * by one withouth having any other boxes behind/above/right of the
 * current box.
 */

void checksol(allinfo *a, box *f, box *l) {
	box *i, *j, *m;
	boolean stillboxes, foundextreme, extreme;
	int b;

	for (i = f, m = l + 1; i != m; i++) {
		if (!i->k)
			continue; /* box currently not chosen */
		for (j = f; j != m; j++) {
			if (i == j)
				continue;
			if (i->no == j->no)
				error("duplicated box %d\n", i->no);
			if (!j->k)
				continue;
			if (i->bno != j->bno)
				continue;
			if ((i->x + i->w > j->x) && (j->x + j->w > i->x) && (i->y + i->h
					> j->y) && (j->y + j->h > i->y) && (i->z + i->d > j->z)
					&& (j->z + j->d > i->z)) {
				error("overlap box %d,%d: [%d,%d,%d] [%d,%d,%d]", i->no, j->no,
						i->w, i->h, i->d, j->w, j->h, j->d);
			}
		}
	}

	if (a->packtype != ROBOT)
		return;

	/* check if robot packable */
	for (b = 1; b <= a->z; b++) {
		for (;;) {
			stillboxes = FALSE;
			foundextreme = TRUE;
			for (i = f; i <= m; i++) {
				if ((i->bno == b) && (i->k == 1))
					stillboxes = TRUE;
			}
			if (!stillboxes)
				break;

			foundextreme = FALSE;
			for (i = f; i <= m; i++) {
				if ((i->bno != b) || (i->k != 1))
					continue;
				extreme = TRUE;
				for (j = f; j <= m; j++) {
					if (j == i)
						continue;
					if ((j->bno != b) || (!j->k != 1))
						continue;
					if ((i->x < j->x + j->w) || (i->y < j->y + j->h) || (i->z
							< j->z + j->d))
						extreme = FALSE;
				}
				if (extreme) {
					i->k = -1;
					foundextreme = TRUE;
				}
			}
			if (!foundextreme)
				break;
		}

		if (!foundextreme) {
			for (i = f; i <= m; i++) {
				if (i->bno == b) {
					printf("no %d (%d %d %d) (%d %d %d) bin %d k %d\n", i->no,
							i->x, i->y, i->z, i->w, i->h, i->d, i->bno, i->k);
				}
			}
			error("not robot packable");
		}
		/* restore k values */
		for (i = f; i <= m; i++)
			if (i->k == -1)
				i->k = 1;
	}
}

/* ======================================================================
 savesol
 ====================================================================== */

/* save an updated solution, checking its validity */

void savesol(allinfo *a, box *f, box *l, ntype z) {
	box *i, *k, *m;

	/* first check validity */
	if (z >= a->z)
		error("not improved");
	for (i = f, m = l + 1; i != m; i++) {
		if ((1 <= i->bno) && (i->bno <= z))
			continue;
		error("illegal bin %d, box %d", i->bno, i->no);
	}

	/* now do the saving */
	a->z = z;
	for (i = f, k = a->fopt, m = l + 1; i != m; i++, k++)
		*k = *i;
	for (i = a->fclosed, m = a->lclosed + 1; i != m; i++, k++)
		*k = *i;
	for (i = a->fopt, m = a->lopt + 1; i != m; i++)
		i->k = TRUE;
	if (DIF(a->fopt,k-1) != a->n)
		error("not correct amount of boxes");
	checksol(a, a->fopt, a->lopt);
}

/* ======================================================================
 isortincr
 ====================================================================== */

/* A specialized routine for sorting integers in increasing order. */
/* qsort could be used equally well, but this routine is faster. */

void isortincr(int *f, int *l) {
	register int mi;
	register int *i, *j, *m;
	register int d;

	d = l - f + 1;
	if (d < 1)
		error("negative interval in isortincr");
	if (d == 1)
		return;
	m = f + d / 2;
	if (*f > *m)
		SWAPINT(f, m);
	if (d > 2) {
		if (*m > *l) {
			SWAPINT(m, l);
			if (*f > *m)
				SWAPINT(f, m);
		}
	}
	if (d <= 3)
		return;
	mi = *m;
	i = f;
	j = l;
	for (;;) {
		do
			i++;
		while (*i < mi);
		do
			j--;
		while (*j > mi);
		if (i > j)
			break;
			else
				SWAPINT(i, j);
	}
	isortincr(f, i - 1);
	isortincr(i, l);
}

/* ======================================================================
 isortdecr
 ====================================================================== */

/* A specialized routine for sorting integers in decreasing order. */
/* qsort could be used equally well, but this routine is faster. */

void isortdecr(int *f, int *l) {
	register int mi;
	register int *i, *j, *m;
	register int d;

	d = l - f + 1;
	if (d < 1)
		error("negative interval in isortdecr");
	if (d == 1)
		return;
	m = f + d / 2;
	if (*f < *m)
		SWAPINT(f, m);
	if (d > 2) {
		if (*m < *l) {
			SWAPINT(m, l);
			if (*f < *m)
				SWAPINT(f, m);
		}
	}
	if (d <= 3)
		return;
	mi = *m;
	i = f;
	j = l;
	for (;;) {
		do
			i++;
		while (*i > mi);
		do
			j--;
		while (*j < mi);
		if (i > j)
			break;
			else
				SWAPINT(i, j);
	}
	isortdecr(f, i - 1);
	isortdecr(i, l);
}

/* ======================================================================
 psortdecr
 ====================================================================== */

/* A specialized routine for sorting extreme points according to decreasing */
/* y-coordinate (decreasing x-coordinate in case of ties) */

void psortdecr(point *f, point *l) {
	register point mi;
	register point *i, *j, *m;
	register int d, r;

	d = l - f + 1;
	if (d <= 1)
		return;
	m = f + d / 2;
	if (DF(*f,*m) < 0)
		SWAPP(f,m);
	if (d == 2)
		return;
	if (DF(*m,*l) < 0) {
		SWAPP(m,l);
		if (DF(*f,*m) < 0)
			SWAPP(f,m);
	}
	if (d <= 3)
		return;
	mi = *m;
	i = f;
	j = l;
	for (;;) {
		do
			i++;
		while (DF(*i,mi) > 0);
		do
			j--;
		while (DF(*j,mi) < 0);
		if (i > j)
			break;
			else
				SWAPP(i, j);
	}
	psortdecr(f, i - 1);
	psortdecr(i, l);
}

/* **********************************************************************
 **********************************************************************
 Lower Bounds
 **********************************************************************
 ********************************************************************** */

/* ======================================================================
 bound_zero
 ====================================================================== */

/* The continuous bound L_0 */

int bound_zero(allinfo *a, box *f, box *l) {
	box *i, *m;
	stype vsum, lb;

	vsum = 0;
	for (i = f, m = l + 1; i != m; i++)
		vsum += i->vol;
	lb = (stype) ceil(vsum / (double) a->BVOL);
	return lb;
}

/* ======================================================================
 rotate_solution
 ====================================================================== */

/* rotates the solution. After 3 rotations we return to original problem */

void rotate_solution(allinfo *a, box *f, box *l) {
	register box *i, *m;
	register itype w, x;

	for (i = f, m = l + 1; i != m; i++) {
		w = i->w;
		i->w = i->h;
		i->h = i->d;
		i->d = w;
		x = i->x;
		i->x = i->y;
		i->y = i->z;
		i->z = x;
	}
}

/* ======================================================================
 rotate_problem
 ====================================================================== */

/* rotates the dimensions by one step */

void rotate_problem(allinfo *a, box *f, box *l) {
	register box *i, *m;
	register itype w, x;

	for (i = f, m = l + 1; i != m; i++) {
		w = i->w;
		i->w = i->h;
		i->h = i->d;
		i->d = w;
		x = i->x;
		i->x = i->y;
		i->y = i->z;
		i->z = x;
	}
	w = a->W;
	a->W = a->H;
	a->H = a->D;
	a->D = w;
}

/* ======================================================================
 choose_boxes
 ====================================================================== */

/* returns a set of boxes with w > W2 and d > D2. This set is used in */
/* bound_one */

void choose_boxes(allinfo *a, box *f, box *l, int W2, int D2, box *fbox,
		box **lbox) {
	box *i, *k, *m;

	for (i = f, m = l + 1, k = fbox; i != m; i++) {
		if ((i->w > W2) && (i->d > D2)) {
			*k = *i;
			k++;
		}
	}
	*lbox = k - 1;
}

/* ======================================================================
 find_plist
 ====================================================================== */

/* returns a zero-terimanted list of distinct dimensions */

void find_plist(box *fbox, box *lbox, itype M, int dim, int *pl) {
	register box *i, *m;
	register int *k, *j, *l;

	i = fbox;
	m = lbox + 1;
	k = pl;
	switch (dim) {
	case WDIM:
		for (; i != m; i++) {
			if (i->w <= M) {
				*k = i->w;
				k++;
			}
		}
		break;
	case HDIM:
		for (; i != m; i++) {
			if (i->h <= M) {
				*k = i->h;
				k++;
			}
		}
		break;
	case DDIM:
		for (; i != m; i++) {
			if (i->d <= M) {
				*k = i->d;
				k++;
			}
		}
		break;
	}
	if (k == pl) {
		*k = 0;
		return;
	}
	isortincr(pl, k - 1); /* sort the dimensions */
	for (j = pl + 1, l = pl; j != k; j++) { /* remove duplicates */
		if (*j != *l) {
			l++;
			*l = *j;
		}
	}
	l++;
	*l = 0;
}

/* ======================================================================
 bound_one
 ====================================================================== */

/* Derive bound L_1 for a fixed dimension */

int bound_one_x(allinfo *a, box *f, box *l) {
	register box *i, *m;
	register itype H, H2, h;
	register int p, j1, j2, j3, j2h, j2hp, j3h;
	int *pp, lb, lb_one, alpha, beta;
	box fbox[MAXBOXES], *lbox;
	int plist[MAXBOXES];

	if (l == f - 1)
		return 0;
	lb = 1;
	H = a->H;
	H2 = H / 2;
	choose_boxes(a, f, l, a->W / 2, a->D / 2, fbox, &lbox);
	if (lbox == fbox - 1) { /* empty */
		return lb;
	}

	find_plist(fbox, lbox, H2, HDIM, plist);
	for (pp = plist; *pp != 0; pp++) {
		p = *pp;
		j1 = j2 = j3 = j2h = j2hp = j3h = 0;
		for (i = fbox, m = lbox + 1; i != m; i++) {
			h = i->h;
			if (h > H - p)
				j1++;
			if ((H - p >= h) && (h > H2)) {
				j2++;
				j2h += h;
				j2hp += (H - h) / p;
			}
			if ((H2 >= h) && (h >= p)) {
				j3++;
				j3h += h;
			}
		}
		alpha = (int) ceil((j3h - (j2 * H - j2h)) / (double) H);
		beta = (int) ceil((j3 - j2hp) / (double) (H / p));
		if (alpha < 0)
			alpha = 0;
		if (beta < 0)
			beta = 0;
		lb_one = j1 + j2 + MAXIMUM(alpha, beta);
		if (lb_one > lb)
			lb = lb_one;
	}
	return lb;
}

/* Derive bound L_1 as the best of all L_1 bounds for three rotations */

int bound_one(allinfo *a, box *f, box *l) {
	int i, lb, lbx;

	lb = 0;
	for (i = WDIM; i <= DDIM; i++) {
		lbx = bound_one_x(a, f, l);
		if (lbx > lb)
			lb = lbx;
		rotate_problem(a, f, l);
	}
	return lb;
}

/* ======================================================================
 bound_two
 ====================================================================== */

/* Derive bound L_2 for a fixed dimension */

int bound_two_x(allinfo *a, box *f, box *l) {
	register box *i, *m;
	register itype W, H, D, w, h, d, W2, D2;
	register int p, q, k1h, k23v;
	int hlb1, lb, lb1, lbx, fract;
	int plist[MAXBOXES], qlist[MAXBOXES];
	int *qq, *pp;
	double WD, BVOL;

	/* derive bound_one */
	lb = lb1 = bound_one_x(a, f, l);
	W = a->W;
	H = a->H;
	D = a->D;
	hlb1 = H * lb1;
	W2 = W / 2;
	D2 = D / 2;
	WD = W * (double) D;
	BVOL = a->BVOL;

	/* run through all values of p, q */
	find_plist(f, l, W2, WDIM, plist);
	find_plist(f, l, D2, DDIM, qlist);
	for (pp = plist; *pp != 0; pp++) {
		p = *pp;
		for (qq = qlist; *qq != 0; qq++) {
			q = *qq;
			k1h = k23v = 0;
			for (i = f, m = l + 1; i != m; i++) {
				w = i->w;
				h = i->h;
				d = i->d;
				if ((w > W - p) && (d > D - q)) {
					k1h += h;
					continue;
				}
				if ((w >= p) && (d >= q)) {
					k23v += i->vol;
				}
			}
			fract = (int) ceil((k23v - (hlb1 - k1h) * WD) / BVOL);
			if (fract < 0)
				fract = 0;
			lbx = lb1 + fract;
			if (lbx > lb)
				lb = lbx;
		}
	}
	return lb;
}

/* Derive bound L_2 as the best of all L_2 bounds for three rotations */

int bound_two(allinfo *a, box *f, box *l) {
	int i, lb, lbx;

	lb = 0;
	for (i = WDIM; i <= DDIM; i++) {
		lbx = bound_two_x(a, f, l);
		if (lbx > lb)
			lb = lbx;
		rotate_problem(a, f, l);
	}
	return lb;
}

/* **********************************************************************
 **********************************************************************
 heuristic filling
 **********************************************************************
 ********************************************************************** */

/* ======================================================================
 onelayer
 ====================================================================== */

/* Fill a layer of depth f->d by arranging the boxes in a number of */
/* vertical shelfs. The boxes $i$ packed are assigned coordinates */
/* (i->x, i->y) and the field i->k is set to the argument d (layer no). */

void onelayer(allinfo *a, box *f, box *m, box *l, int d) {
	int s, t;
	itype r; /* remaining width */
	itype width[MAXBOXES], height[MAXBOXES], x[MAXBOXES];
	box *i;

	qsort(f, DIF(f,m), sizeof(box), (funcptr) hcomp);
	r = a->W;
	x[0] = 0;
	width[0] = 0;
	height[0] = 0;
	for (s = 1, i = f; i != l + 1; s++) {
		x[s] = x[s - 1] + width[s - 1];
		height[s] = 0;
		width[s] = i->w;
		if (width[s] > r)
			width[s] = r;
		r -= width[s];
		for (; i != l + 1; i++) {
			for (t = s; t != 0; t--) {
				if (i->w <= width[t]) {
					if (height[t] + i->h <= a->H) {
						i->y = height[t];
						i->x = x[t];
						i->k = d;
						height[t] += i->h;
						break;
					}
				}
			}
			if ((t == 0) && (r > 0))
				break; /* new strip */
		}
	}
}

/* ======================================================================
 countarea
 ====================================================================== */

/* Select a subset of the boxes such that the selected boxes have a total */
/* area of two times the face of a bin (the parameter: barea) */

box *countarea(box *f, box *l, stype barea) {
	box *i;
	stype area, d;

	for (area = 0, i = f; i != l + 1; i++) {
		d = i->h * (ptype) i->w;
		area += d;
		if (area > 2 * barea)
			return i - 1;
	}
	return l;
}

/* ======================================================================
 remboxes
 ====================================================================== */

/* Remove the boxes which were chosen for a layer, i.e., where i->k != 0. */
/* The depth of the layer is set equal to the deepest box chosen. */

box *remboxes(box *f, box *l, itype *depth) {
	box *i, *j;
	itype d;

	for (i = f, j = l, d = 0; i != j + 1;) {
		if (i->k) {
			if (i->d > d)
				d = i->d;
			i++;
		} else {
			SWAP(i,j);
			j--;
		}
	}
	*depth = d;
	return i;
}

/* ======================================================================
 assignboxes
 ====================================================================== */

/* Assign z-coordinates to the boxes, once they layers have been combined */
/* to individual bins by solving a 1-dimensional Bin-packing Problem. */

void assignboxes(heurpair *t, heurpair *u, ntype maxbin, box *f, box *l) {
	box *i, *m;
	heurpair *h;
	itype b, z;

	/* derive z-coordinates for each layer */
	for (b = 1; b <= maxbin; b++) {
		z = 0;
		for (h = t; h <= u; h++)
			if (h->bno == b) {
				h->z = z;
				z += h->d;
			}
	}

	for (i = f, m = l + 1; i != m; i++) {
		h = t + i->k - 1;
		i->z = h->z;
		i->bno = h->bno;
	}
}

/* ======================================================================
 onedim_binpack
 ====================================================================== */

/* One-dimensional bin-packing algorithm. In each iteration, the next */
/* box is assigned to every open bin as well as to a new bin. The  */
/* algorithm terminates when MAXBPP iterations have been performed, */
/* returning the heuristic solution found. */

void onedim_binpack(heurpair *i, heurpair *f, heurpair *l, int *b, int bno,
		itype *z) {
	int j, *bc;
	heurpair *k;

	bpiterat++;
	if (bpiterat > MAXBPP)
		return;
	if (bno >= *z)
		return; /* no hope of improvement */

	if (i > l) {
		*z = bno;
		for (k = f; k <= l; k++)
			k->bno = k->b;
	} else {
		for (j = 0; j < bno; j++) {
			bc = b + j;
			if (i->d <= *bc) {
				*bc -= i->d;
				i->b = j + 1;
				onedim_binpack(i + 1, f, l, b, bno, z);
				*bc += i->d;
			}
		}
		bc = b + bno;
		*bc -= i->d;
		i->b = bno + 1;
		onedim_binpack(i + 1, f, l, b, bno + 1, z);
		*bc += i->d;
	}
}

/* ======================================================================
 dfirst_heuristic
 ====================================================================== */

/* Heuristic algorithm for the 3D BPP. A number of layers are constructed */
/* using the shelf approach to pack every layer. Then the individual layers */
/* are combined to bins by solving a 1D BPP defined in the layer depths. */

void dfirst_heuristic(allinfo *a) {
	box *j, *f, *l, *m;
	int i, n, h, b[MAXBOXES];
	heurpair t[MAXBOXES];
	itype d, z;

	/* initialize boxes */
	for (j = a->fbox, m = a->lbox + 1; j != m; j++) {
		j->bno = j->x = j->y = j->z = 0;
		j->k = FALSE;
	}

	/* fill layer one by one */
	for (f = a->fbox, l = a->lbox, h = 0;; h++) {
		n = DIF(f,l);
		if (n == 0)
			break;
		qsort(f, n, sizeof(box), (funcptr) dcomp);
		m = countarea(f, l, a->W * (ptype) a->H);
		onelayer(a, f, m, l, h + 1);
		f = remboxes(f, l, &d);
		t[h].d = d;
		t[h].bno = h + 1; /* initially put layers into separate bins */
		t[h].z = 0;
		t[h].lno = h + 1; /* this ensures fes. solution if terminate */
	}

	/* split into bins by solving 1-dim binpacking */
	for (i = 0; i < h; i++)
		b[i] = a->D; /* all bins are empty */
	qsort(t, h, sizeof(heurpair), (funcptr) xcomp);
	z = h + 1;
	bpiterat = 0;
	onedim_binpack(t, t, t + h - 1, b, 0, &z);
	qsort(t, h, sizeof(heurpair), (funcptr) lcomp); /* order according to lno */

	/* now assign bin number to each boxes */
	assignboxes(t, t + h - 1, z, a->fbox, a->lbox);
	if (z < a->zlayer)
		a->zlayer = z;
	if (a->zlayer < a->z)
		savesol(a, a->fbox, a->lbox, a->zlayer);
}

/* ======================================================================
 dfirst3_heuristic
 ====================================================================== */

/* Call the heuristic dfirst_heuristic, for three different rotations */
/* of the problem */

void dfirst3_heuristic(allinfo *a) {
	int i;
	double t1, t2;

	timer(&t1);
	a->zlayer = a->n; /* very bad incumbent solution */
	for (i = WDIM; i <= DDIM; i++) {
		dfirst_heuristic(a);
		rotate_solution(a, a->fopt, a->lopt);
		rotate_problem(a, a->fbox, a->lbox);
	}
	timer(&t2);
	a->lhtime = t2 - t1;
}

/* **********************************************************************
 **********************************************************************
 fill one 3D bin using GENERAL packing
 **********************************************************************
 ********************************************************************** */

/* ======================================================================
 modifyandpush
 ====================================================================== */

/* Push the relation "rel" between box i and box j to a stack. If "dom"
 * is true, the relation is removed from the domain. If "dom" is false,
 * the relation "rel" is imposed between boxes "i" and "j".
 */

void modifyandpush(int i, int j, int rel, boolean dom) {
	dompos->i = i;
	dompos->j = j;
	dompos->domain = dom;
	if (dom) {
		dompos->relation = rel;
		domain[i][j][rel] = FALSE;
	} else {
		dompos->relation = relation[i][j];
		relation[i][j] = rel;
	}
	dompos++;
	if (dompos == domend)
		error("stack filled\n");
}

/* ======================================================================
 popdomains
 ====================================================================== */

/* Pop all relations between boxes from the stack. The stack is emptied
 * downto the depth given by "pos".
 */

void popdomains(domainpair *pos) {
	for (; dompos != pos;) {
		dompos--;
		if (dompos->domain) {
			domain[dompos->i][dompos->j][dompos->relation] = TRUE;
		} else {
			relation[dompos->i][dompos->j] = dompos->relation;
		}
	}
}

/* ======================================================================
 findcoordinates
 ====================================================================== */

/* Find coordinates of boxes according to the currently pending relations.
 * In principle this can be done by topologically sorting the boxes according
 * to e.g. the left-right relations, and then assigning coordinates from
 * the left-most box to the right-most box.
 *   The following implementation is a simplified version, which runs through
 * all pairs of boxes, and checks whether they satisfy the relation, otherwise
 * moving one of the boxes right, up or behind, according to the given
 * relation. The process is repeated until no pairs of boxes violate a
 * relation. Computational experiments have shown that the present approach
 * for the considered instances is faster than a topological sorting followed
 * by a critical path calculation.
 *   If a box during the process gets moved outsides of the bin, then
 * the algorithm terminates with FALSE. Otherwise TRUE is returned, saying
 * that a feasible packing exists which respects the current relations.
 */

boolean findcoordinates(allinfo *a, int n, box *f) {
	register box *g, *h;
	register int sum;
	int i, j, k, W, H, D;
	boolean changed;
	char *dom, *relij;
	domset *domij;

	/* check if feasible, i.e., at least one choice for each relation */
	W = a->W;
	H = a->H;
	D = a->D;
	for (i = 0; i < n; i++) {
		j = i + 1;
		relij = &(relation[i][j]);
		domij = &(domain[i][j]);
		for (; j < n; j++, relij++, domij++) {
			if (*relij != UNDEF)
				continue;
			dom = *domij;
			if (*dom)
				continue;
			dom++;
			if (*dom)
				continue;
			dom++;
			if (*dom)
				continue;
			dom++;
			if (*dom)
				continue;
			dom++;
			if (*dom)
				continue;
			dom++;
			if (*dom)
				continue;
			return FALSE;
		}
	}

	/* initialize coordinates */
	for (i = 0; i < n; i++) {
		g = f + i;
		g->x = g->y = g->z = 0;
	}

	/* now determine the coordinates */
	a->exacttopo++;
	for (k = 0; k < n; k++) {
		a->exacttopn++;
		changed = FALSE;
		for (i = 0; i < n; i++) {
			g = f + i;
			j = i + 1;
			relij = &(relation[i][j]);
			for (; j < n; j++, relij++) {
				h = f + j;
				switch (*relij) {
				case UNDEF:
					/* do nothing */
					break;
				case LEFT:
					sum = g->x + g->w;
					if (h->x < sum) {
						h->x = sum;
						changed = TRUE;
						if (sum + h->w > W)
							return FALSE;
					}
					break;
				case RIGHT:
					sum = h->x + h->w;
					if (g->x < sum) {
						g->x = sum;
						changed = TRUE;
						if (sum + g->w > W)
							return FALSE;
					}
					break;
				case UNDER:
					sum = g->y + g->h;
					if (h->y < sum) {
						h->y = sum;
						changed = TRUE;
						if (sum + h->h > H)
							return FALSE;
					}
					break;
				case ABOVE:
					sum = h->y + h->h;
					if (g->y < sum) {
						g->y = sum;
						changed = TRUE;
						if (sum + g->h > H)
							return FALSE;
					}
					break;
				case FRONT:
					sum = g->z + g->d;
					if (h->z < sum) {
						h->z = sum;
						changed = TRUE;
						if (sum + h->d > D)
							return FALSE;
					}
					break;
				case BEHIND:
					sum = h->z + h->d;
					if (g->z < sum) {
						g->z = sum;
						changed = TRUE;
						if (sum + g->d > D)
							return FALSE;
					}
					break;
				}
			}
		}
		if (!changed) {
			return TRUE;
		}
	}
	/* there must be a loop in the graph */
	return FALSE;
}

/* ======================================================================
 checkdomain
 ====================================================================== */

/* Temporarily impose the relation "value" between boxes "i" and "j",
 * and check whether a feasible assignment of coordinates exists which
 * respects all currently imposed relations. 
 *   If the relation cannot be satisfied, it is removed from the domain
 * and pushed to a stack, so that it can be restored upon backtracking.
 */

void checkdomain(allinfo *a, int i, int j, int n, box *f, int value) {
	if (domain[i][j][value] == FALSE)
		return; /* not allowed in any case */
	relation[i][j] = value;
	if (findcoordinates(a, n, f) == FALSE) {
		modifyandpush(i, j, value, TRUE);
	}
}

/* ======================================================================
 reducedomain
 ====================================================================== */

/* Constraint propagation algorithm. For each relation in the domain of
 * boxes "i" and "j", check if the relation has the posibility of being
 * satisfied. If some of the relations cannot be satisfied any more, they
 * are removed from the domain (and pushed to a stack, so that they can
 * be restored when the master search algorithm backtracks). If only one
 * relation remains in the domain, the relation is imposed at this node
 * and all descendant nodes.
 */

boolean reducedomain(allinfo *a, int n, box *f) {
	register int i, j, k, l, m;

	m = 0;
	for (i = 0; i < n - 1; i++) {
		for (j = i + 1; j < n - 1; j++) {
			if (relation[i][j] == UNDEF) {
				checkdomain(a, i, j, n, f, LEFT);
				checkdomain(a, i, j, n, f, RIGHT);
				checkdomain(a, i, j, n, f, UNDER);
				checkdomain(a, i, j, n, f, ABOVE);
				checkdomain(a, i, j, n, f, FRONT);
				checkdomain(a, i, j, n, f, BEHIND);
				relation[i][j] = UNDEF;
				for (k = LEFT, l = 0; k < UNDEF; k++) {
					if (domain[i][j][k]) {
						l++;
						m = k;
					}
				}
				if (l == 0)
					return FALSE;
				if (l == 1) {
					modifyandpush(i, j, m, FALSE);
				}
			}
		}
	}
	return TRUE;
}

/* ======================================================================
 recpack
 ====================================================================== */

/* Recursive algorithm based on constraint programming used for assigning
 * relative positions to each pair of boxes. Each pair of boxes initially 
 * has an associated relation with domain LEFT, RIGHT, UNDER, ABOVE, FRONT, 
 * BEHIND. In each iteration of the algorithm a pair of boxes "i" and "j"
 * is assigned the relation "rel". Constraint propagation is then used to 
 * decrease the domains of remaining relations. 
 *   If it is not possible to assign coordinates to the boxes such that the 
 * currently imposed relations between pairs of boxes are respected, we 
 * backtrack.
 *   If each pair of boxes has been assigned a relation, and it is possible
 * to assign coordinates to the boxes such that the currently imposed 
 * relations between pairs of boxes are respected, we save the solution 
 * and return.
 *   Otherwise, constraint propagation is used to decrease the domains
 * of relations corresponding to each pairs of boxes. If a domain only
 * contains a single relation, the relation is fixed.
 *   The recursive step selects the next pair of boxes following "i" and "j"
 * and repeatedly assigns each relation from the domain to the relation 
 * variable.
 */

void recpack(allinfo *a, int i, int j, int n, box *f, int rel) {
	int i1, j1;
	domainpair *pos;
	boolean feas;

	if (stopped)
		return;
	a->iter3d++;
	if ((a->iter3d == a->maxiter) && (a->maxiter != 0))
		terminate = TRUE;
	a->subiterat++;
	if (a->subiterat == IUNIT) {
		a->subiterat = 0;
		a->iterat++;
		check_iterlimit(a->iterat, a->iterlimit);
		check_timelimit(a->timelimit);
	}
	if (terminate)
		return;

	relation[i][j] = rel;

	for (i1 = 0, j1 = 0; i1 != i && j1 != j;) {
		i1++;
		if (i1 >= j1) {
			i1 = 0;
			j1++;
		}
		if (relation[i1][j1] == UNDEF)
			error("relation error %d %d\n", i1, j1);
	}

	feas = findcoordinates(a, n, f);
	if (!feas)
		return;

	if ((i == n - 2) && (j == n - 1)) {
		feasible = TRUE;
		terminate = TRUE;
		memcpy(a->fsol, f, sizeof(box) * n);
		return;
	}

	pos = dompos;
	feas = reducedomain(a, n, f);
	if (feas) {
		i++;
		if (i >= j) {
			i = 0;
			j++;
		}
		bblevel++;
		rel = relation[i][j];
		if (domain[i][j][LEFT])
			recpack(a, i, j, n, f, LEFT);
		if (domain[i][j][RIGHT])
			recpack(a, i, j, n, f, RIGHT);
		if (domain[i][j][UNDER])
			recpack(a, i, j, n, f, UNDER);
		if (domain[i][j][ABOVE])
			recpack(a, i, j, n, f, ABOVE);
		if (domain[i][j][FRONT])
			recpack(a, i, j, n, f, FRONT);
		if (domain[i][j][BEHIND])
			recpack(a, i, j, n, f, BEHIND);
		relation[i][j] = rel;
		bblevel--;
	}
	popdomains(pos);
}

/* ======================================================================
 general_pack
 ====================================================================== */

/* General packing procedure, which tests whether boxes f..l can be packed
 * into a single bin. The algorithm is based on constraint programming, where
 * each pair of boxes initially has an associated relation with domain LEFT, 
 * RIGHT, UNDER, ABOVE, FRONT, BEHIND. Then the recursive algorithm "recpack" 
 * is called, which repeatedly tries to assign the relation a value, using
 * constraint propagation to decrease the domains of remaining boxes.
 */

boolean general_pack(allinfo *a, box *f, box *l) {
	register int i, j, k, n;

	dompos = domstack;
	domend = domstack + STACKDEPTH;
	feasible = FALSE;
	terminate = FALSE;
	bblevel = 1;
	n = l - f + 1;
	if (n > a->exactn)
		a->exactn = n;

	for (i = 0; i < n; i++) {
		for (j = 0; j < n; j++) {
			relation[i][j] = UNDEF;
			for (k = LEFT; k < UNDEF; k++) {
				domain[i][j][k] = TRUE;
			}
		}
	}
	domain[0][1][RIGHT] = FALSE;
	domain[0][1][ABOVE] = FALSE;
	domain[0][1][BEHIND] = FALSE;

	recpack(a, 0, 0, n, f, UNDEF);
	return feasible;
}

/* ======================================================================
 onebin_general
 ====================================================================== */

/* Check if boxes f..l can be packed into a single bin using general 
 * packing. If "fast" is TRUE, the problem is only solved heuristically,
 * hence if the algorithm returns FALSE we cannot rule out the posibility
 * of packing all boxes into the bin.
 */

boolean onebin_general(allinfo *a, box *f, box *l, boolean fast) {
	boolean solution;
	double t1, t2;

	/* check time limit */
	if (stopped)
		return FALSE;

	a->iter3d = 0;
	a->maxiter = (fast ? MAXITER : 0); /* limited or infinitly many */
	a->exactcall++;

	/* calling general pack */
	timer(&t1);
	solution = general_pack(a, f, l);
	if (solution)
		checksol(a, f, l);
	timer(&t2);
	a->genertime += t2 - t1;
	return solution;
}

/* **********************************************************************
 **********************************************************************
 fill one 3D bin using ROBOT packing
 **********************************************************************
 ********************************************************************** */

/* ======================================================================
 envelope
 ====================================================================== */

/* Find the two-dimensional envelope of the boxes given by extreme */
/* points f to l. */

void envelope(point *f, point *l, point *s1, point **sm, itype W, itype H,
		itype D, itype RW, itype RH, itype cz, point **ll, int *nz, stype *area) {
	register point *i, *s, *t;
	register itype x, xx, y, z, ix, iy, iz, mz;
	register stype sum;

	/* find corner points and area */
	x = xx = z = 0;
	y = H;
	sum = 0;
	mz = D;
	for (i = t = f, s = s1; i != l; i++) {
		iz = i->z;
		if (iz <= cz)
			continue;
		if (iz < mz)
			mz = iz; /* find minimum next z coordinate */
		ix = i->x;
		if (ix <= x) {
			if (iz > z) {
				*t = *i;
				t++;
			}
			continue;
		}
		iy = i->y;
		if ((x <= RW) && (iy <= RH)) {
			s->x = x;
			s->y = iy;
			s->z = cz;
			s++;
			sum += (x - xx) * (ptype) y;
			y = iy;
			xx = x;
		}
		x = ix;
		z = iz;
		*t = *i;
		t++;
	}
	if (y != 0)
		sum += (W - xx) * (ptype) y;
	*sm = s - 1;
	*area = sum;
	*nz = mz;
	*ll = t - 1;
}

/* ======================================================================
 checkdom
 ====================================================================== */

/* The 3D envelope is found by deriving a number of 2D envelopes. This */
/* may however introduce some "false" corner points, which are marked */
/* by the following algorithm. */

void checkdom(point *s1, point *sl, point *sm) {
	register point *s, *t, *u;

	if (sl == s1 - 1)
		return;
	for (s = s1, t = sl + 1, u = sm + 1; t != u; t++) {
		while (s->x < t->x) {
			s++;
			if (s > sl)
				return;
		}
		if ((s->x == t->x) && (s->y == t->y))
			t->z = 0;
	}
}

/* ======================================================================
 removedom
 ====================================================================== */

/* Remove "false" corner points marked by algorithm "checkdom". */

point *removedom(point *s1, point *sl) {
	register point *i, *m, *k;

	for (i = k = s1, m = sl + 1; i != m; i++) {
		if (i->z == 0)
			continue;
		*k = *i;
		k++;
	}
	return k - 1;
}

/* ======================================================================
 initboxes
 ====================================================================== */

/* Initialize boxes. Already placed boxes define extreme points, and thus */
/* they are placed into the list fc,..,lc. Not placed boxes are used to */
/* derive minimum dimensions. */

void initboxes(box *f, box *l, point *fc, point **lc, int *minw, int *minh,
		int *mind) {
	register box *j, *m;
	register point *k;
	register int minx, miny, minz;

	minx = *minw;
	miny = *minh;
	minz = *mind;
	for (j = f, k = fc, m = l + 1; j != m; j++) {
		if (j->k) { /* defines an extreme box */
			k->x = j->x + j->w;
			k->y = j->y + j->h;
			k->z = j->z + j->d;
			k++;
		} else { /* free box */
			if (j->w < minx)
				minx = j->w;
			if (j->h < miny)
				miny = j->h;
			if (j->d < minz)
				minz = j->d;
		}
	}
	*minw = minx;
	*minh = miny;
	*mind = minz;
	*lc = k - 1;
}

/* ======================================================================
 findplaces
 ====================================================================== */

/* Find all corner points, where a new box may be placed as well as the */
/* volume of the "envelope" occupied by already placed boxes. Already */
/* placed boxes are given by f,..,l, while a list of possible placings */
/* is returned in s1,..,sm. The return value of the procedure is an upper */
/* bound on the possible filling of this bin. */

stype findplaces(allinfo *a, box *f, box *l, point *s1, point **sm, stype fill) {
	register point *k;
	int minw, minh, mind, W, H, D, RW, RH, z, zn;
	point *sk, *lc, *sl, *st, *s0;
	stype vol, area;
	point fc[MAXBOXES + 1];

	/* select boxes which are chosen, and find min dimensions of unchosen */
	minw = W = a->W;
	minh = H = a->H;
	mind = D = a->D;
	initboxes(f, l, fc, &lc, &minw, &minh, &mind);

	/* sort the boxes according to max y (first) max x (second) */
	if (lc >= fc)
		psortdecr(fc, lc); /* order decreasing */

	/* for each z-coordinate find the 2D envelope */
	vol = 0;
	sl = s1 - 1;
	sk = s1;
	s0 = NULL;
	RW = W - minw;
	RH = H - minh;
	lc++;
	k = lc;
	k->x = W + 1;
	k->y = 0;
	k->z = a->D + 1;
	for (z = 0; z != D; z = zn) {
		/* find 2D envelope for all boxes which cover *z */
		envelope(fc, lc + 1, sl + 1, &st, W, H, D, RW, RH, z, &lc, &zn, &area);
		if (zn + mind > D)
			zn = D; /* nothing fits between zn and D */
		vol += area * (ptype) (zn - z); /* update volume */
		checkdom(sk, sl, st); /* check for dominance */
		sk = sl + 1;
		sl = st;
		if (z == 0)
			s0 = sl;
	}
	*sm = removedom(s0 + 1, sl);
	return fill + (a->BVOL - vol); /* bound is curr filling + all free vol */
}

/* ======================================================================
 branch
 ====================================================================== */

/* Recursive algorithm for solving a knapsack filling of a single bin. */
/* In each iteration, the set of feasible positions for placing a new */
/* box is found, and an upper bound on the filling is derived. If the */
/* bound indicates that an improved solution still may be obtained, every */
/* box is tried to be placed at every feasible position, before calling */
/* the algorithm recursively. */

void branch(allinfo *a, box *f, box *l, int miss, stype fill) {
	register box *i;
	register point *s;
	int d;
	stype bound;
	point *sl, s1[9 * MAXBOXES];

	if (stopped)
		return;
	a->iter3d++;
	if ((a->iter3d == a->maxiter) && (a->maxiter != 0))
		terminate = TRUE;
	if (a->iter3d % 1000 == 0)
		check_timelimit(a->timelimit);
	a->subiterat++;
	if (a->subiterat == IUNIT) {
		a->subiterat = 0;
		a->iterat++;
		check_iterlimit(a->iterat, a->iterlimit);
	}
	if (terminate)
		return;

	/* find min/max dimensions of remaining boxes */
	if (miss == 0) {
		/* none left -> good: save solution */
		memcpy(a->fsol, f, sizeof(box) * DIF(f,l));
		a->maxfill = a->BVOL;
		terminate = TRUE;
		a->miss = miss;
	} else {
		/* check if better filling */
		if (fill > a->maxfill) {
			memcpy(a->fsol, f, sizeof(box) * DIF(f,l));
			a->maxfill = fill;
			a->miss = miss;
		}

		/* find bound and positions to place new boxes */
		bound = findplaces(a, f, l, s1, &sl, fill);

		if (bound > a->maxfill) {
			/* for each position in S, try to place an box there */
			for (s = s1; s != sl + 1; s++) {
				for (i = f, d = 0; i != l + 1; i++) {
					if (i->k)
						continue; /* already chosen */

					/* see if box fits at position s */
					if ((int) (s->x) + i->w > a->W)
						continue;
					if ((int) (s->y) + i->h > a->H)
						continue;
					if ((int) (s->z) + i->d > a->D)
						continue;

					/* place box and call recursively */
					i->k = TRUE;
					i->x = s->x;
					i->y = s->y;
					i->z = s->z;
					branch(a, f, l, miss - 1, fill + i->vol);
					i->k = FALSE;
					i->x = i->y = i->z = 0;
					d++;
					if (d == a->mcut)
						break; /* terminate after mcut branches */
					if (terminate)
						break;
				}
			}
		}
	}
}

/* ======================================================================
 onebin_robot
 ====================================================================== */

/* Knapsack filling of a single bin. The following procedure initializes */
/* some variables before calling the recursive branch-and-bound algorithm. */

boolean onebin_robot(allinfo *a, box *f, box *l, boolean fast) {
	box *i, *j, *m;
	stype vol;
	double t1, t2;
	int n;

	/* initialize boxes */
	for (i = f, m = l + 1, vol = 0; i != m; i++) {
		i->x = 0;
		i->y = 0;
		i->z = 0;
		i->k = 0;
		vol += i->vol;
	}

	/* try to fill one bin with boxes [f,l] */
	n = DIF(f,l);
	a->iter3d = 0;
	a->maxfill = vol - 1; /* lower bound equal to all minus one */
	a->miss = n;
	a->maxiter = (fast ? MAXITER : 0); /* limited or infinitly many */
	a->mcut = 0; /* try all branches */
	terminate = FALSE;

	/* calling branch-and-bound algorithm */
	timer(&t1);
	branch(a, f, l, n, 0);
	timer(&t2);
	a->robottime += t2 - t1;

	/* copy solution */
	if (a->maxfill == a->BVOL) { /* NB: maxfill = BVOL, when optimal found */
		for (i = a->fsol, j = f, m = l + 1; j != m; i++, j++)
			*j = *i;
		return TRUE;
	}
	return FALSE;
}

/* ======================================================================
 mcut_heuristic
 ====================================================================== */

/* Knapsack filling of a single bin, solved heuristically. The heuristic */
/* is based on the exact algorithm for knapsack filling a single bin, where */
/* only a limited number of sub-nodes are considered at every branching node */
/* (the so-called m-cut approach) */

void mcut_heuristic(allinfo *a) {
	box *f, *l, *i, *j, *m;
	int b, n;

	/* initialize boxes */
	for (i = a->fbox, m = a->lbox + 1; i != m; i++) {
		i->bno = i->x = i->y = i->z = 0;
		i->k = FALSE;
	}

	/* fill bins one by one */
	f = a->fbox;
	l = a->lbox;
	for (b = 1;; b++) {
		/* fill one bin */
		for (i = f; i <= l; i++)
			i->k = FALSE; /* box not chosen */
		a->iter3d = 0;
		a->maxfill = 0;
		a->miss = DIF(f,l);
		a->maxiter = 5 * MAXITER;
		terminate = FALSE;
		n = DIF(f,l);
		a->mcut = 2;
		if (n < 15)
			a->mcut = 3;
		if (n < 10)
			a->mcut = 4;
		branch(a, f, l, n, 0);

		/* copy solution */
		for (i = a->fsol, j = f, m = l + 1; j != m; i++, j++)
			*j = *i;

		/* remove chosen boxes */
		for (i = f; i <= l;)
			if (i->k) {
				i->bno = b;
				SWAP(i,l);
				l--;
			} else
				i++;

		/* check if finished */
		if (l == f - 1)
			break;
	}
	if (b < a->zmcut)
		a->zmcut = b; /* save solution */
	if (a->zmcut < a->z)
		savesol(a, a->fbox, a->lbox, a->zmcut);
}

/* ======================================================================
 mcut3_heuristic
 ====================================================================== */

/* Knapsack filling of a single bin, solved heuristically. Three different */
/* rotations of the problem are considered, and the best found solution */
/* is selected. */

void mcut3_heuristic(allinfo *a) {
	int i;
	double t1, t2;

	timer(&t1);
	a->zmcut = a->n; /* very bad lower bound */
	for (i = WDIM; i <= DDIM; i++) {
		mcut_heuristic(a);
		rotate_solution(a, a->fopt, a->lopt);
		rotate_problem(a, a->fbox, a->lbox);
	}
	timer(&t2);
	a->mhtime = t2 - t1;
}

/* **********************************************************************
 **********************************************************************
 branch-and-bound for 3D bin-packing problem
 **********************************************************************
 ********************************************************************** */

/* ======================================================================
 fits
 ====================================================================== */

/* The routine "fitsm" checks whether a given subset of boxes fits into */
/* a single bin. To improve performance, specialized algorithms are derived */
/* for cases with two to three boxes */

boolean fits2(box *i, box *j, itype W, itype H, itype D) {
	/* all coordinates are initialized to zero, so just adjust changes! */
	/* the 2-box solution is always guillotine cuttable */
	if (i->w + j->w <= W) {
		j->x = i->w;
		return TRUE;
	}
	if (i->h + j->h <= H) {
		j->y = i->h;
		return TRUE;
	}
	if (i->d + j->d <= D) {
		j->z = i->d;
		return TRUE;
	}
	return FALSE;
}

boolean fits2p(box *i, box *j, itype W, itype H, itype D) {
	if (i->w + j->w <= W)
		return TRUE;
	if (i->h + j->h <= H)
		return TRUE;
	if (i->d + j->d <= D)
		return TRUE;
	return FALSE;
}

boolean fits3(box *i, box *j, box *k, itype W, itype H, itype D, int packtype) {
	box *t;
	itype w, h, d, r;

	/* all coordinates are initialized to zero, so just adjust changes! */
	/* the 3-box solution can either be cut by guillotine cuts */
	for (r = 1; r <= 3; r++) {
		/* cut (i,j) and (k) in one of three dimensions */
		w = W - k->w;
		h = H - k->h;
		d = D - k->d;
		if ((i->w <= w) && (j->w <= w) && fits2(i, j, w, H, D)) {
			k->x = w;
			return TRUE;
		}
		if ((i->h <= h) && (j->h <= h) && fits2(i, j, W, h, D)) {
			k->y = h;
			return TRUE;
		}
		if ((i->d <= d) && (j->d <= d) && fits2(i, j, W, H, d)) {
			k->z = d;
			return TRUE;
		}
		t = i;
		i = j;
		j = k;
		k = t;
	}
	/* or by a sticky arrangement which is not guillotine packable */
	if ((packtype != GENERAL) && (packtype != ROBOT))
		return FALSE;

	/* (xi,yi,zi) = (0,0,0); (xj,yj,zj) = (wi,0,0); (xk,yk,zk) = (0,hi,dj) */
	if ((i->w + j->w <= W) && (i->h + k->h <= H) && (j->d + k->d <= D)) {
		j->x = i->w;
		k->y = i->h;
		k->z = j->d;
		return TRUE;
	}
	/* (xi,yi,zi) = (0,0,0); (xj,yj,zj) = (wk,0,di); (xk,yk,zk) = (0,hi,0) */
	if ((j->w + k->w <= W) && (i->h + k->h <= H) && (i->d + j->d <= D)) {
		j->x = k->w;
		j->z = i->d;
		k->y = i->h;
		return TRUE;
	}
	/* (xi,yi,zi) = (0,0,0); (xj,yj,zj) = (0,hi,dk); (xk,yk,zk) = (wi,0,0) */
	if ((i->w + k->w <= W) && (i->h + j->h <= H) && (k->d + j->d <= D)) {
		j->y = i->h;
		j->z = k->d;
		k->x = i->w;
		return TRUE;
	}
	/* (xi,yi,zi) = (0,0,0); (xj,yj,zj) = (0,hi,0); (xk,yk,zk) = (wj,0,di) */
	if ((j->w + k->w <= W) && (i->h + j->h <= H) && (k->d + i->d <= D)) {
		j->y = i->h;
		k->x = j->w;
		k->z = i->d;
		return TRUE;
	}
	/* (xi,yi,zi) = (0,0,0); (xj,yj,zj) = (wi,0,0); (xk,yk,zk) = (0,hj,di) */
	if ((i->w + j->w <= W) && (j->h + k->h <= H) && (i->d + k->d <= D)) {
		j->x = i->w;
		k->y = j->h;
		k->z = i->d;
		return TRUE;
	}
	/* (xi,yi,zi) = (0,0,0); (xj,yj,zj) = (0,0,di); (xk,yk,zk) = (wi,hj,0) */
	if ((i->w + k->w <= W) && (j->h + k->h <= H) && (i->d + j->d <= D)) {
		j->z = i->d;
		k->x = i->w;
		k->y = j->h;
		return TRUE;
	}
	return FALSE;
}

boolean fitsm(allinfo *a, box *t, box *k, boolean fast) {
	boolean fits;
	ntype lb;

	lb = bound_two(a, t, k);
	if (lb > 1)
		return FALSE;
	a->exfill++;
	fits = FALSE;
	if (a->packtype == GENERAL)
		fits = onebin_general(a, t, k, fast);
	if (a->packtype == ROBOT)
		fits = onebin_robot(a, t, k, fast);
	return fits;
}

/* ======================================================================
 onebin_decision
 ====================================================================== */

/* The following procedure checks whether a new box "j" fits into the the */
/* bin "bno" (together with already placed boxes in the bin). If the answer */
/* is "no" then definitly it is not possible to place the box into the bin. */

boolean onebin_decision(allinfo *a, box *j, int bno) {
	register box *i, *k, *m;
	box t[MAXBOXES];
	boolean fits;
	int size;

	for (i = a->fbox, m = j, k = t - 1; i != m; i++) {
		if (i->bno == bno) {
			k++;
			*k = *i;
			k->x = k->y = k->z = 0;
			k->ref = i;
		}
	}
	k++;
	*k = *j;
	k->x = k->y = k->z = 0;
	k->ref = j;
	k->k = TRUE;

	size = DIF(t,k);
	switch (size) {
	case 0:
		error("no boxes in onebin_decision");
	case 1:
		fits = TRUE;
		k->x = k->y = k->z = 0;
		break;
	case 2:
		fits = fits2(t, k, a->W, a->H, a->D);
		break;
	case 3:
		fits = fits3(t, t + 1, k, a->W, a->H, a->D, a->packtype);
		break;
	default:
		fits = fitsm(a, t, k, FALSE);
		break;
	}
	if (size <= 3) {
		a->subiterat++;
		if (a->subiterat == IUNIT) {
			a->subiterat = 0;
			a->iterat++;
			check_iterlimit(a->iterat, a->iterlimit);
		}
	}

	if (fits) {
		/* copy solution back */
		for (i = t, m = k + 1; i != m; i++) {
			k = i->ref;
			k->x = i->x;
			k->y = i->y;
			k->z = i->z;
			k->k = TRUE;
		}
	}
	return fits;
}

/* ======================================================================
 onebin_heuristic
 ====================================================================== */

/* This is a heuristic version of "onebin_decision". If the answer is "yes" */
/* then a heuristic solution has been found where boxes f,..,l fit into a  */
/* bin. If the answer is "no" then a filling may still be possible, but it */
/* was not found by the heuristic. */

boolean onebin_heuristic(allinfo *a, box *f, box *l) {
	box *i, *m;
	boolean fits;

	for (i = f, m = l + 1; i != m; i++) {
		i->x = i->y = i->z = 0;
	}
	switch (DIF(f,l)) {
	case 0:
		error("no boxes in onebin_heuristic");
	case 1:
		fits = TRUE;
		break;
	case 2:
		fits = fits2(f, l, a->W, a->H, a->D);
		break;
	case 3:
		fits = fits3(f, f + 1, l, a->W, a->H, a->D, a->packtype);
		break;
	default:
		fits = fitsm(a, f, l, TRUE);
		break;
	}
	return fits;
}

/* ======================================================================
 try_close
 ====================================================================== */

/* A bin may be closed if no more boxes fit into it. The present version */
/* uses a more advanced criterion: First, the set of boxes which fit into */
/* bin "bno" is derived. This is done, by testing whether each box (alone) */
/* fits together with the already placed boxes in the bin. Having derived */
/* the set of additional boxes that (individually) fits into the bin, it is */
/* tested whether a solution exists where all the additional boxes are */
/* placed in the bin. If this is the case, we have found a optimal placing */
/* of the additional boxes, and thus we may close the bin. */

boolean try_close(allinfo *a, box **curr, ntype bno, box *oldf, box **oldl,
		box **oldlc, ntype *oldnoc, boolean *oldclosed, int level) {
	register box *j, *m, *k, *r, *i, *s;
	register stype vol;
	box f[MAXBOXES];
	ntype h, n, b;
	boolean didclose, fits;

	if (level > MAXCLOSE)
		return FALSE;
	i = *curr;
	didclose = FALSE;
	for (b = 1; b <= bno; b++) {
		if (i > a->lbox)
			break;
		if (a->closed[b])
			continue;
		for (j = a->fbox, m = i, k = f, n = 0, vol = 0; j != m; j++) {
			if (j->bno == b) {
				*k = *j;
				k->ref = j;
				k++;
				n++;
				vol += j->vol;
			}
		}
		if (n == 0)
			error("bin with no boxes");
		if (vol < a->BVOL / 2)
			continue;
		for (j = i, h = 0, m = a->lbox + 1; j != m; j++) {
			if ((j->no > a->n) || (j->no < 1))
				error("bad no");
			fits = onebin_decision(a, j, b);
			if (fits) {
				*k = *j;
				k->ref = j;
				k++;
				h++;
				vol += j->vol;
			}
			if (vol > a->BVOL)
				break;
		}
		if (vol > a->BVOL)
			continue;
		if (onebin_heuristic(a, f, k - 1)) {
			if (!didclose) { /* take backup of table when first bin closed */
				memcpy(oldclosed, a->closed, sizeof(boolean) * (bno + 1));
				memcpy(oldf, a->fbox, sizeof(box) * DIF(a->fbox,a->lbox));
				*oldl = a->lbox;
				*oldlc = a->lclosed;
				*oldnoc = a->noc;
			}
			a->closed[b] = TRUE;
			s = a->lclosed;
			didclose = TRUE;
			a->noc++;
			if (a->noc > a->maxclose)
				a->maxclose = a->noc;
			for (j = f; j != k; j++) {
				r = j->ref;
				r->bno = b;
				r->k = TRUE;
				r->x = j->x;
				r->y = j->y;
				r->z = j->z;
			}
			for (j = k = a->fbox, m = a->lbox + 1; j != m; j++) {
				if (j->bno == b) {
					s++;
					*s = *j;
				} else {
					*k = *j;
					k++;
				}
			}
			a->lbox = k - 1;
			a->lclosed = s;
			i -= n; /* reposition current box */
		}
	}
	*curr = i;
	return didclose;
}

/* ======================================================================
 free_close
 ====================================================================== */

/* Reopen a closed bin when backtracking. */

void free_close(allinfo *a, ntype bno, box *oldf, box *oldl, box *oldlc,
		ntype oldnoc, boolean *oldclosed) {
	a->lbox = oldl;
	a->lclosed = oldlc;
	a->noc = oldnoc;
	memcpy(a->fbox, oldf, sizeof(box) * DIF(a->fbox,oldl));
	memcpy(a->closed, oldclosed, sizeof(boolean) * (bno + 1));
}

/* ======================================================================
 rec_binpack
 ====================================================================== */

/* Recursive algorithm for 3D Bin-packing Problem. In each iteration, the */
/* next box "i" is assigned to every open bin, as well as to a new bin. */

void rec_binpack(allinfo *a, box *i, int bno, ntype lb, int level) {
	box of[MAXBOXES], *ol, *ox;
	boolean ocl[MAXBOXES];
	ntype b, oc;
	boolean more;

	if (bno >= a->z)
		return; /* used too many bins */
	if (a->z == a->lb)
		return; /* optimal solution found */
	a->subnodes++;
	if (a->subnodes == IUNIT) {
		a->subnodes = 0;
		a->nodes++;
	}
	check_nodelimit(a->nodes, a->nodelimit);
	check_iterlimit(a->iterat, a->iterlimit);
	check_timelimit(a->timelimit);
	if (stopped)
		return;

	if (i == a->lbox + 1) {
		/* all boxes assigned, must be better solution */
		savesol(a, a->fbox, a->lbox, bno);
	} else {
		more = try_close(a, &i, bno, of, &ol, &ox, &oc, ocl, level);
		if (i == a->lbox + 1) { /* all boxes went into closed bins */
			savesol(a, a->fbox, a->lbox, bno);
		} else {
			if (more)
				lb = a->noc + bound_two(a, a->fbox, a->lbox);
			if (lb < a->z) {
				for (b = 1; b <= bno; b++) {
					if (a->closed[b])
						continue; /* cannot add to closed bin */
					if (onebin_decision(a, i, b)) {
						i->bno = b;
						rec_binpack(a, i + 1, bno, lb, level + 1);
						i->bno = 0;
					}
				}
				i->bno = bno + 1;
				i->x = i->y = i->z = 0;
				a->closed[i->bno] = FALSE;
				rec_binpack(a, i + 1, bno + 1, lb, level + 1);
				i->bno = 0;
			}
		}
		/* restore */
		if (more)
			free_close(a, bno, of, ol, ox, oc, ocl);
	}
}

/* **********************************************************************
 **********************************************************************
 Main procedure
 **********************************************************************
 ********************************************************************** */

/* ======================================================================
 clearboxes
 ====================================================================== */

void clearboxes(allinfo *a) {
	box *i, *m;

	for (i = a->fbox, m = a->lbox + 1; i != m; i++) {
		i->x = i->y = i->z = i->bno = 0;
		i->k = FALSE;
		i->vol = VOL(i);
	}
	/* sort nonincreasing volume */
	qsort(a->fbox, (m - a->fbox), sizeof(box), (funcptr) vcomp);
}

/* ======================================================================
 copyboxes
 ====================================================================== */

void copyboxes(allinfo *a, int *w, int *h, int *d, int W, int H, int D) {
	box *i, *m;
	int k;

	for (i = a->fbox, m = a->lbox + 1, k = 0; i != m; i++, k++) {
		i->no = k + 1;
		i->w = w[k];
		i->h = h[k];
		i->d = d[k];
		if ((w[k] < 1) || (w[k] > W))
			error("bad w\n");
		if ((h[k] < 1) || (h[k] > H))
			error("bad h\n");
		if ((d[k] < 1) || (d[k] > D))
			error("bad d\n");
	}

	clearboxes(a);
}

/* ======================================================================
 returnboxes
 ====================================================================== */

void returnboxes(allinfo *a, int *x, int *y, int *z, int *bno) {
	box *i, *m;
	int k;

	for (i = a->fopt, m = a->lopt + 1; i != m; i++) {
		k = i->no - 1;
		x[k] = i->x;
		y[k] = i->y;
		z[k] = i->z;
		bno[k] = i->bno;
	}
}

/* ======================================================================
 binpack3d
 ====================================================================== */

void binpack3d(int n, int W, int H, int D, int *w, int *h, int *d, int *x,
		int *y, int *z, int *bno, int *lb, int *ub, int nodelimit,
		int iterlimit, int timelimit, int *nodeused, int *iterused,
		int *timeused, int packingtype) {
	allinfo a;
	box t0[MAXBOXES], t1[MAXBOXES], t2[MAXBOXES], t3[MAXBOXES];
	boolean cl[MAXBOXES];

	/* start the timer */
	timer(NULL);
	stopped = FALSE;

	/* copy info to a structure */
	if ((packingtype != ROBOT) && (packingtype != GENERAL))
		error("bad packtype");
	if (n + 1 > MAXBOXES)
		error("too big instance");
	a.n = n;
	a.W = W;
	a.H = H;
	a.D = D;
	a.fbox = t0;
	a.lbox = a.fbox + a.n - 1;
	a.fsol = t1;
	a.lsol = a.fsol + a.n - 1;
	a.fopt = t2;
	a.lopt = a.fopt + a.n - 1;
	a.fclosed = t3;
	a.lclosed = a.fclosed - 1;
	a.noc = 0;
	a.closed = cl;
	a.BVOL = W * (ptype) H * D;
	a.packtype = packingtype;
	a.maxfill = 0;
	a.exfill = 0;
	a.nodelimit = 0;
	a.iterlimit = 0;
	a.timelimit = 0;
	a.nodes = 0;
	a.subnodes = 0;
	a.iterat = 0;
	a.subiterat = 0;
	a.didpush = 0;
	a.maxclose = 0;
	a.genertime = 0;
	a.robottime = 0;
	a.z = a.n + 1;

	/* copy boxes to internal structure */
	copyboxes(&a, w, h, d, W, H, D);

	/* find bounds */
	a.bound0 = bound_zero(&a, a.fbox, a.lbox);
	a.bound1 = bound_one(&a, a.fbox, a.lbox);
	a.bound2 = bound_two(&a, a.fbox, a.lbox);
	a.lb = a.bound2;

	/* find heuristic solution */
	dfirst3_heuristic(&a);

	/* initialize search limits for exact search */
	a.nodelimit = nodelimit;
	a.iterlimit = iterlimit;
	a.timelimit = timelimit;

	/* outer tree enummeration */
	clearboxes(&a); /* clear positions */
	rec_binpack(&a, a.fbox, 0, a.lb, 1);
	timer(&(a.time));

	/* check found solution */
	/* checksol(&a, a.fopt, a.lopt); */

	/* copy boxes back to arrays */
	returnboxes(&a, x, y, z, bno);
	*ub = a.z;
	*lb = (stopped ? a.lb : a.z);
	*nodeused = a.nodes;
	*iterused = a.iterat;
	*timeused = a.time * 1000;
}