diff --git a/src/set/ibex_NodeType.cpp b/src/set/ibex_NodeType.cpp
index d379266a0..24918e1b0 100644
--- a/src/set/ibex_NodeType.cpp
+++ b/src/set/ibex_NodeType.cpp
@@ -20,7 +20,7 @@
 // ========================================
 
 namespace ibex {
-
+// relic from irregular set, need for SetBisect constructor
 NodeType operator|(NodeType x, NodeType y) {
 	switch (x) {
 	case IN : {
@@ -78,6 +78,219 @@ NodeType operator|(NodeType x, NodeType y) {
 	}
 }
 
+
+
+// define operator or between two nodetype, assuming the two boxes associated are equal
+NodeType uni(NodeType x,NodeType y) {
+	switch (x) {
+	case IN : {
+		return IN;
+		
+		}
+	break;
+	case OUT : {
+		return y;
+		}
+	break;
+	case UNK : {
+		switch(y) {
+		case IN :     return IN;
+		case UNK_IN:  return UNK_IN;
+		case OUT :    return UNK;
+		case UNK :	  return UNK;
+		case UNK_OUT: return UNK_OUT;
+		default:      return UNK_IN_OUT;
+		}
+	}
+	break;
+	case UNK_IN : {
+		switch(y) {
+		case IN :
+		case UNK_IN:
+		case UNK :    return UNK_IN;
+		default:      return UNK_IN_OUT;
+		}
+	}
+	break;
+	case UNK_OUT : {
+		switch(y) {
+		case OUT :
+		case UNK_OUT:
+		case UNK :    return UNK_OUT;
+		default:      return UNK_IN_OUT;
+		}
+	}
+	break;
+	default :
+		//Warning: IN_TMP not considered here.
+		return UNK_IN_OUT;
+	}
+
+}
+
+// define operator or between two nodetype, assuming box of nodetype y is subset of box of nodetype x;
+NodeType uni_in(NodeType x, NodeType y) {
+	switch (x) {
+	case IN : {
+		switch(y) {
+		case IN :        return IN;
+		case OUT :		 return IN;
+		case UNK :		 return IN;
+		case UNK_OUT:
+		case UNK_IN_OUT: return UNK_IN_OUT;
+		default :        return UNK_IN;
+		}
+	}
+	break;
+	case OUT : {
+		switch(y) {
+		case OUT :       return OUT;
+		case IN :		 return UNK;
+		case UNK :		 return UNK;
+		case UNK_IN:
+		case UNK_IN_OUT: return UNK_IN_OUT;
+		default :        return UNK_OUT;
+		}
+	}
+	break;
+	case UNK : {
+		switch(y) {
+		return UNK;
+		}
+	}
+	break;
+	/*case UNK_IN : {
+		switch(y) {
+		case IN :
+		case UNK_IN:
+		case UNK :    return UNK_IN;
+		default:      return UNK_IN_OUT;
+		}
+	}
+	break;
+	case UNK_OUT : {
+		switch(y) {
+		case OUT :
+		case UNK_OUT:
+		case UNK :    return UNK_OUT;
+		default:      return UNK_IN_OUT;
+		}
+	}
+	break;*/
+	default :
+		//Warning: IN_TMP not considered here.
+		return UNK_IN_OUT;
+		
+	}
+}
+
+// define operator and between two nodetype, assuming the two boxes associated are equal
+NodeType inte(NodeType x, NodeType y) {
+	switch (x) {
+	case IN : {
+		return y;
+		}
+	break;
+	case OUT : {
+		return OUT;
+		}
+	break;
+	case UNK : {
+		switch(y) {
+		case IN :     return UNK;
+		case UNK_IN:  return UNK_IN;
+		case OUT :    return OUT;
+		case UNK :	  return UNK;
+		case UNK_OUT: return UNK_OUT;
+		default:      return UNK_IN_OUT;
+		}
+	}
+	break;
+	case UNK_IN : {
+		switch(y) {
+		case IN :
+		case UNK_IN:
+		case UNK :    return UNK_IN;
+		default:      return UNK_IN_OUT;
+		}
+	}
+	break;
+	case UNK_OUT : {
+		switch(y) {
+		case OUT :
+		case UNK_OUT:
+		case UNK :    return UNK_OUT;
+		default:      return UNK_IN_OUT;
+		}
+	}
+	break;
+	default :
+		//Warning: IN_TMP not considered here.
+		return UNK_IN_OUT;
+	}
+}
+
+
+// define operator and between two nodetype, assuming box of nodetype y is subset of box of nodetype x;
+NodeType inte_in(NodeType x, NodeType y) {
+	switch (x) {
+	case IN : {
+		switch(y) {
+		case IN :        return IN;
+		case OUT :		 return UNK;
+		case UNK :		 return UNK;
+		case UNK_OUT:
+		case UNK_IN_OUT: return UNK_IN_OUT;
+		default :        return UNK_IN;
+		}
+	}
+	break;
+	case OUT : {
+		switch(y) {
+		case OUT :       return OUT;
+		case IN :		 return OUT;
+		case UNK :		 return OUT;
+		case UNK_IN:
+		case UNK_IN_OUT: return UNK_IN_OUT;
+		default :        return UNK_OUT;
+		}
+	}
+	break;
+	case UNK : {
+		switch(y) {
+		case IN :	  return UNK;
+		case UNK_IN:  return UNK_IN;
+		case OUT :	  return UNK;
+		case UNK_OUT: return UNK_OUT;
+		case UNK:     return UNK;
+		default:      return UNK_IN_OUT;
+		}
+	}
+	break;
+	case UNK_IN : {
+		switch(y) {
+		case IN :
+		case UNK_IN:
+		case UNK :    return UNK_IN;
+		default:      return UNK_IN_OUT;
+		}
+	}
+	break;
+	case UNK_OUT : {
+		switch(y) {
+		case OUT :
+		case UNK_OUT:
+		case UNK :    return UNK_OUT;
+		default:      return UNK_IN_OUT;
+		}
+	}
+	break;
+	default :
+		//Warning: IN_TMP not considered here.
+		return UNK_IN_OUT;
+	}
+}
+
 bool certainly_contains_in(NodeType x) {
 	return x==IN || x==UNK_IN || x==UNK_IN_OUT;
 }
diff --git a/src/set/ibex_NodeType.h b/src/set/ibex_NodeType.h
index 8d060ede1..7456e64f4 100644
--- a/src/set/ibex_NodeType.h
+++ b/src/set/ibex_NodeType.h
@@ -44,9 +44,24 @@ typedef enum { __IBEX_IN__,
  * \see SetBisect constructor.
  */
 NodeType operator|(NodeType x, NodeType y);
-
 /**
- * \brief False only if the subtree contains no inner point
+ * \brief return value of union of two equal boxes
+ */
+NodeType uni(NodeType x,NodeType y);
+/**
+ * \brief return value of union of two boxes, assuming the on associated to y is include in the one associated to x
+ */
+NodeType uni_in(NodeType x, NodeType y);
+/**
+ * \brief return value of inter of two equal boxes
+ */
+NodeType inte(NodeType x, NodeType y);
+/**
+ * \brief return value of inter of two boxes, assuming the on associated to y is include in the one associated to x
+ */
+NodeType inte_in(NodeType x, NodeType y);
+/**
+ * \brief return value of inter of two boxes, assuming the on associated to y is include in the one associated to x
  */
 bool possibly_contains_in(NodeType x);
 
diff --git a/src/set/ibex_SetBisect.cpp b/src/set/ibex_SetBisect.cpp
index ead14b71d..858147302 100644
--- a/src/set/ibex_SetBisect.cpp
+++ b/src/set/ibex_SetBisect.cpp
@@ -26,14 +26,15 @@ using namespace std;
 
 namespace ibex {
 
-SetBisect::SetBisect(int var, double pt, SetNode* left, SetNode* right) : SetNode(left->status | right->status), var(var), pt(pt), left(left), right(right) {
-	// a bisectNode with two subnodes of same status IN or OUT should not exist
-	// (automatically compacted as a leaf node) but two subnodes IN_TMP can be
-	// created by a leaf with IN_TMP status (when it auto-splits in inter function)
-	assert(left->status>=UNK || left->status!=right->status);
+SetBisect::SetBisect(int var, double pt, SetNode* left, SetNode* right) : SetNode(UNK_IN_OUT,NULL), var(var), pt(pt), left(left), right(right) {
+	
 }
 
-SetBisect::SetBisect(int var, double pt) : SetNode(UNK), var(var), pt(pt), left(NULL), right(NULL) {
+SetBisect::SetBisect(int var, double pt, SetNode* left, SetNode* right,SetNode* father) : SetNode(UNK_IN_OUT,father), var(var), pt(pt), left(left), right(right) {
+	
+}
+
+SetBisect::SetBisect(int var, double pt) : SetNode(UNK_IN_OUT), var(var), pt(pt), left(NULL), right(NULL) {
 
 }
 
@@ -43,87 +44,271 @@ SetBisect::~SetBisect() {
 }
 
 bool SetBisect::is_leaf() const {
-	return false;
+    return false;
+}
+
+SetNode * SetBisect::copy() const {
+	SetBisect * n = new SetBisect(this->var,this->pt,this->left->copy(),this->right->copy());
+	n->left->father = n;
+	n->right->father = n;
+	return n;
+}
+
+
+void SetBisect::inter(NodeType x) {
+
+	left->inter(x); // apply inter to left and right until reach the leaves
+	right->inter(x);
 }
 
-SetNode* SetBisect::sync(const IntervalVector& nodebox, const IntervalVector& x, NodeType x_status, double eps) {
-	assert(x_status<=UNK);
+void SetBisect::_union(NodeType x) {
+	left->_union(x);
+	right->_union(x);
+}
 
-	if (x_status==UNK) {
-		return this;
+void SetBisect::oper(SetNode * node,bool op) {
+	if(node->is_leaf()) // apply status to leaves of the subtree which have this as root
+	{
+		if(op && node->status != IN) 
+		{
+			left->inter(node->status);
+			right->inter(node->status);
+		}
+		else if(!op && node->status!= OUT)
+		{
+			left->_union(node->status);
+			right->_union(node->status);
+		}
 	}
-	else if (nodebox.is_subset(x)) {
-		if (x_status==IN && !possibly_contains_in(status)) throw NoSet();
-		if (x_status==OUT && !possibly_contains_out(status)) throw NoSet();
-
-		delete this; // warning: suicide
-		return new SetLeaf(x_status);
-	} else {
-		left = left->sync(left_box(nodebox), x, x_status, eps);
-		right = right->sync(right_box(nodebox), x, x_status, eps);
-		// status of children may have changed --> try merge
-		return try_merge();
+	else
+	{
+		SetBisect * other = (SetBisect*) node;  // node is assume to be either a leaf or a bisect, able to cast node in bisect as it ain't a leaf
+		left->oper(other->left,op);
+		right->oper(other->right,op);
 	}
 }
 
-SetNode* SetBisect::sync_rec(const IntervalVector& nodebox, Sep& sep, double eps) {
-	left = left->sync(left_box(nodebox), sep, eps);
-	right = right->sync(right_box(nodebox), sep, eps);
-	// status of children may have changed --> try merge
-	return try_merge();
+void SetBisect::operator_ir(const IntervalVector& box,const IntervalVector& subbox, NodeType val,bool op, double eps)
+{
+
+	if(box.is_subset(subbox)) // compute intersection of leaves of this root with val
+	{
+		if(op && val != IN)
+			this->inter(val);
+		else if(!op && val != OUT)
+			this->_union(val);
+	}
+	else
+	{
+		if(this->left_box(box).overlaps(subbox)) //  need to apply function on left only if leftbox ovelaps subbox
+			left->operator_ir(this->left_box(box),subbox,val,op,eps);
+		if(this->right_box(box).overlaps(subbox)) //  need to apply function on right only if rightbox overlaps subbox
+			right->operator_ir(this->right_box(box),subbox,val,op,eps);
+	}
 }
 
+void SetBisect::operator_ir(const IntervalVector& box,const IntervalVector& subbox, NodeType valin,NodeType valout,bool op, double eps)
+{
 
-SetNode* SetBisect::inter(const IntervalVector& nodebox, const IntervalVector& x, NodeType x_status, double eps) {
-	assert(x_status<=UNK);
+	if(box.is_subset(subbox)) // apply valin on the leaves of this root
+	{
+		if(op && valin != IN)
+			this->inter(valin);
+		else if(!op && valin != OUT)
+			this->_union(valin);
+	}
+	else if(box.overlaps(subbox))
+	{
+		if(this->left_box(box).overlaps(subbox))
+			left->operator_ir(this->left_box(box),subbox,valin,valout,op,eps);
+		else {								// leftbox is disjoint of subbox, apply valout on the leaves of this->left root
+			if(op && valout != IN)
+				left->inter(valout);
+			else if (!op && valout != OUT)
+				left->_union(valout);}
+		if(this->right_box(box).overlaps(subbox))
+			right->operator_ir(this->right_box(box),subbox,valin,valout,op,eps);
+		else{								// rightbox is disjoint of subbox, apply valout on the leaves of this->right root
+			if(op && valout != IN)
+				right->inter(valout);
+			else if (!op && valout != OUT)
+				right->_union(valout);}
+	}
+	else{   								// this is disjoint of subbox, apply valout on the leaves of this root
+		if(op && valout != IN)
+			this->inter(valout);
+		else if (!op && valout != OUT)
+			this->_union(valout);}
+}
 
-	// no: keep subnodes informed that x_status is IN (their status can changed from IN_TMP to IN)
-	//	if (x_status==IN) {
-	//		return this;
-	//	}
+// Uncomment this function and comment the one above to use fakeBranch method
+/*void SetBisect::operator_ir(const IntervalVector& box,const IntervalVector& subbox, NodeType val,bool op, double eps)
+{
 
-	// in comment because certainly_contains_in does not take into account IN_TMP
-//	if (x_status==UNK && !certainly_contains_in(status)) {
-//		return this;
-//	}
+	if(left == NULL)
+	{
+		assert(right == NULL);
+		left = this->fakeLeaf(left_box(box),subbox,val,op,eps);
+		left->father = this;
+		right = this->fakeLeaf(right_box(box),subbox,val,op,eps);
+		right->father = this;
+		status = UNK_IN_OUT;
+	}
+	else
+	{
+		if(box.is_subset(subbox))
+		{
+			if(op)
+				this->inter(val);
+			else
+				this->_union(val);
+		}
+		else
+		{
+			if(this->left_box(box).overlaps(subbox))
+				left->operator_ir(this->left_box(box),subbox,val,op,eps);
+			else
+				return;
+			if(this->right_box(box).overlaps(subbox))
+				right->operator_ir(this->right_box(box),subbox,val,op,eps);
+			else
+				return;
+		}
+	}
+}*/
 
-	// certainly_contains_out in comment because does not take into account IN_TMP
-    if ((x_status==OUT /*|| !certainly_contains_out(status)*/) && nodebox.is_subset(x)) {
-		delete this; // warning: suicide
-		return new SetLeaf(x_status); // either OUT or UNK
-	} else {
-		left = left->inter(left_box(nodebox), x, x_status, eps);
-		right = right->inter(right_box(nodebox), x, x_status, eps);
-		// status of children may have changed --> try merge
-		return try_merge();
+SetNode * SetBisect::fakeLeaf(const IntervalVector& box,const IntervalVector& subbox,NodeType val,bool op, double eps) {
+    //assert(left == NULL && right == NULL);
+	if(box.is_subset(subbox))
+		if(op)
+			return new SetLeaf(inte(status,val));
+		else
+			return new SetLeaf(uni(status,val));
+	else if(box.overlaps(subbox))
+	{
+		if(box.max_diam()<=eps)
+			return new SetLeaf(inte_in(status,val));
+		else
+		{
+			int var = box.extr_diam_index(false); 
+			double pt = box[var].mid();
+			SetBisect * node = new SetBisect(var,pt,NULL,NULL);
+			node->status  = status;
+			node->left = node->fakeLeaf(node->left_box(box),subbox,val,op,eps);
+			node->right = node->fakeLeaf(node->right_box(box),subbox,val,op,eps);
+			node->left->father = node;
+			node->right->father = node;
+			node->status = UNK_IN_OUT;
+			return node;
+		}
 	}
+	else
+		return new SetLeaf(status);
 }
+void SetBisect::cleave(const IntervalVector& box, Sep& sep, const double eps) {
 
-SetNode* SetBisect::inter_rec(const IntervalVector& nodebox, Sep& sep, double eps) {
-	left = left->inter(left_box(nodebox), sep, eps);
-	right = right->inter(right_box(nodebox), sep, eps);
-	// status of children may have changed --> try merge
-	return try_merge();
-}
+	/*IntervalVector box1(box);
+	IntervalVector box2(box);
 
-SetNode* SetBisect::union_(const IntervalVector& nodebox, const IntervalVector& x, NodeType x_status, double eps) {
-	assert(x_status<=UNK);
+	sep.separate(box1,box2);
+	if(box1.is_empty())
+		this->inter(OUT);
+	else if(box2.is_empty())
+		this->_union(IN);
+	else // continu until box1 and box2 are disjoint
+	{
+		left->cleave(left_box(box),sep,eps);
+		right->cleave(right_box(box),sep,eps);
+	}*/
+	
+	status = UNK_IN_OUT;
+	IntervalVector box1(box);
+	IntervalVector box2(box);
 
-	if (x_status>IN) {
-		return this;
+	sep.separate(box1,box2);
+	if(box1.is_empty()) // all the box can be set to OUT
+	{
+		this->inter(OUT);
 	}
+	else if(box2.is_empty())// all the box can be set to IN
+	{
+		this->_union(IN);
+	}
+	else
+	{
 
-    if (nodebox.is_subset(x)) {
-		delete this; // warning: suicide
-		return new SetLeaf(IN);
-	} else {
-		left = left->union_(left_box(nodebox), x, x_status, eps);
-		right = right->union_(right_box(nodebox), x, x_status, eps);
-		// status of children may have changed --> try merge
-		return try_merge();
+		IntervalVector* restout;
+    	int n=0;
+    	if(box1!=box)
+    		box.diff(box1,restout);
+    	for(int i = 0;i<n;i++)
+    		this->operator_ir(box,restout[i],OUT,true,eps); // add in box
+    	IntervalVector* restin;
+    	n = 0;
+    	if(box2!= box)
+    		n=box.diff(box2,restin);
+    	for(int i = 0;i<n;i++)
+    		this->operator_ir(box,restin[i],IN,false,eps); // add out box
+    	this->gatherTo(false,this);
+    	if(!left->is_leaf() || left->status == UNK )//&& left_box(box).max_diam()>eps)
+    		left->cleave(left_box(box),sep,eps);
+    	if(!right->is_leaf() || right->status == UNK )//&& right_box(box).max_diam()>eps)
+    		right->cleave(right_box(box),sep,eps);
 	}
+	
+}
+void SetBisect::gather(bool go_up) {
+		if(((left->is_leaf()&&right->is_leaf()) || go_up)&&left->status==right->status) // leaves are reached, or climbing up the tree toward the root
+		{
+			status = left->status;  // this bisect get a leaf status IN OUT or UNK
+			if(father != NULL) // check if father is not the root of SetInterval
+				father->gather(true);
+		}
+		else if(!go_up) // go down the tree to reach the leaves
+		{
+			left->gather(false);
+			right->gather(false);
+		}
 }
 
+void SetBisect::gatherTo(bool go_up, SetNode * branch) {
+		if(((left->is_leaf()&&right->is_leaf()) || go_up)&&left->status==right->status) // leaves are reached, or climbing up the tree toward the root
+		{
+			status = left->status;  // this bisect get a leaf status IN OUT or UNK
+			if(father != branch) // check if father is not the root of SetInterval
+				father->gather(true);
+		}
+		else if(!go_up) // go down the tree to reach the leaves
+		{
+			left->gather(false);
+			right->gather(false);
+		}
+}
+
+void SetBisect::cutDeadBranch() {
+	if(left->status ==right->status && left->status < 3) // if status of left and right are equal and has got a leaf status, meaning all leaves have the same status
+	{
+		SetBisect* bfather = (SetBisect*) father;
+		if(bfather->left == this)
+			bfather->left = new SetLeaf(status,father); // this is now a leaf
+		else
+			bfather->right = new SetLeaf(status,father); // this is now a leaf
+		delete this;  // delete former this which was a bisect
+	}
+	else // go down the tree until reach leaves or bisect that can be set a leaf
+	{
+		left->cutDeadBranch();
+		right->cutDeadBranch();
+	}
+}
+
+
+void SetBisect::checkFat() {
+	if(father == NULL)
+		cout<<"issue, invalid father"<<endl;
+	left->checkFat();
+	right->checkFat();
+}
 void SetBisect::visit_leaves(leaf_func func, const IntervalVector& nodebox) const {
 	left->visit_leaves(func, left_box(nodebox));
 	right->visit_leaves(func, right_box(nodebox));
@@ -136,15 +321,13 @@ void SetBisect::print(ostream& os, const IntervalVector& nodebox, int shift) con
 	right->print(os, right_box(nodebox), shift+2);
 }
 
-void SetBisect::set_in_tmp() {
-	left->set_in_tmp();
-	right->set_in_tmp();
+void SetBisect::setFathers() {
+	left->father = this;
+	right->father = this;
+	left->setFathers();
+	right->setFathers();
 }
 
-void SetBisect::unset_in_tmp() {
-	left->unset_in_tmp();
-	right->unset_in_tmp();
-}
 
 IntervalVector SetBisect::left_box(const IntervalVector& nodebox) const {
 	IntervalVector leftbox(nodebox);
@@ -160,14 +343,6 @@ IntervalVector SetBisect::right_box(const IntervalVector& nodebox) const {
 	return rightbox;
 }
 
-SetNode* SetBisect::try_merge() {
-	// the case left=right=UNK may happen.
-	if (left->status<=UNK && left->status==right->status) {
-		NodeType s=left->status;
-		delete this;
-		return new SetLeaf(s);
-	} else
-		return this;
-}
 
 } // namespace ibex
+
diff --git a/src/set/ibex_SetBisect.h b/src/set/ibex_SetBisect.h
index e9b9eadd3..fb0073e38 100644
--- a/src/set/ibex_SetBisect.h
+++ b/src/set/ibex_SetBisect.h
@@ -30,6 +30,8 @@ class SetBisect : public SetNode {
 	 */
 	SetBisect(int var, double pt, SetNode* left, SetNode* right);
 
+	SetBisect(int var, double pt, SetNode* left, SetNode* right,SetNode* father);
+
 	/**
 	 * \brief Delete this
 	 */
@@ -39,19 +41,37 @@ class SetBisect : public SetNode {
 	virtual bool is_leaf() const;
 
 	/** \see SetNode */
-	virtual SetNode* sync(const IntervalVector& nodebox, const IntervalVector& x, NodeType x_status, double eps);
+	virtual SetNode * copy() const;
+
+	/** \see SetNode */
+	virtual void inter(NodeType x_status);
+
+	/** \see SetNode */
+	virtual void _union(NodeType x);
+
+	/** \see SetNode */
+	virtual void oper(SetNode * other,bool op);
+
+	/** \see SetNode */
+	virtual void operator_ir(const IntervalVector& box,const IntervalVector& nodebox, NodeType val,bool op, double eps);
+
+	/** \see SetNode */
+	virtual void operator_ir(const IntervalVector& box,const IntervalVector& subbox, NodeType valin,NodeType valout, bool op,double eps);
 
 	/** \see SetNode */
-	virtual SetNode* inter(const IntervalVector& nodebox, const IntervalVector& x, NodeType x_status, double eps);
+	virtual SetNode * fakeLeaf(const IntervalVector& box,const IntervalVector& subbox,NodeType val,bool op, double eps);
 
 	/** \see SetNode */
-	virtual SetNode* sync_rec(const IntervalVector& nodebox, Sep& sep, double eps);
+	virtual void cleave(const IntervalVector& box, Sep& sep, double eps);
 
 	/** \see SetNode */
-	virtual SetNode* inter_rec(const IntervalVector& nodebox, Sep& sep, double eps);
+	virtual void gather(bool go_up);
 
 	/** \see SetNode */
-	virtual SetNode* union_(const IntervalVector& nodebox, const IntervalVector& x, NodeType x_status, double eps);
+	virtual void cutDeadBranch();
+
+	/** \see SetNode */
+	virtual void checkFat();
 
 	/** \see SetNode */
 	virtual void visit_leaves(leaf_func func, const IntervalVector& nodebox) const;
@@ -60,13 +80,20 @@ class SetBisect : public SetNode {
 	virtual void print(std::ostream& os, const IntervalVector& nodebox, int shift) const;
 
 	/** \see SetNode */
-	virtual void set_in_tmp();
+	virtual void setFathers();
+
+
 
-	/** \see SetNode */
-	virtual void unset_in_tmp();
 
+
+	/**
+	 * \brief Rerturn left_box according to var, pt and nodebox
+	 */
 	IntervalVector left_box(const IntervalVector& nodebox) const;
 
+	/**
+	 * \brief Rerturn left_box according to var, pt and nodebox
+	 */
 	IntervalVector right_box(const IntervalVector& nodebox) const;
 
 //protected:
@@ -76,7 +103,6 @@ class SetBisect : public SetNode {
 	// partial initialization used by SetInterval::load only
 	SetBisect(int var, double pt);
 
-	SetNode* try_merge();
 
 	int var;
 	double pt;
diff --git a/src/set/ibex_SetInterval.cpp b/src/set/ibex_SetInterval.cpp
index 0dc1c1d5f..691bc13b2 100644
--- a/src/set/ibex_SetInterval.cpp
+++ b/src/set/ibex_SetInterval.cpp
@@ -27,38 +27,93 @@ SetInterval::SetInterval(const char* filename) : root(NULL), eps(-1), bounding_b
 	load(filename);
 }
 
+SetInterval::SetInterval(const SetInterval& set) : eps(set.eps), bounding_box(set.bounding_box)
+{
+	root = set.root->copy();
+	root->father = NULL;
+}
+
+
 bool SetInterval::is_empty() const {
 	return root==NULL;
 }
 
-void SetInterval::sync(Sep& sep) {
-	try {
-		root = root->sync(bounding_box, sep, eps);
-	} catch(NoSet& e) {
-		delete root;
-		root = NULL;
-		throw e;
-	}
-}
+void SetInterval::cutRoot()
+{
+	NodeType stat = root->status;
+	delete root;
+	int var = bounding_box.extr_diam_index(false);
+	double pt = bounding_box[var].mid();
+	SetBisect * broot = new SetBisect(var,pt,new SetLeaf(stat),new SetLeaf(stat));
+	broot->left->father = broot;
+	broot->right->father = broot;
+	root = broot;
 
+}
 void SetInterval::contract(Sep& sep) {
-	root->set_in_tmp();
-	root = root->inter(bounding_box, sep, eps);
-	root->unset_in_tmp();
+	// case root is leaf create issue for SetLeaf recursive function that need the father
+	// root is redefined as SetBisect pointing on two leaves with value of the root
+	if(root->is_leaf())
+		this->cutRoot();
+	root->cleave(bounding_box, sep,eps);
+	this->gather();
+	root->father = NULL;
 }
 
 SetInterval& SetInterval::operator&=(const SetInterval& set) {
-	root->set_in_tmp();
-	root = root->inter(bounding_box, set.root, set.bounding_box, eps);
-	root->unset_in_tmp();
+	assert(bounding_box== set.bounding_box); // root must be the same as set interval are regular
+	// case root is leaf create issue for SetLeaf recursive function that need the father
+	// root is redefined as SetBisect pointing on two leaves with value of the root
+	if(root->is_leaf())
+		this->cutRoot();
+	root->oper(set.root,true);
+	this->gather();
 	return *this;
 }
 
+SetInterval& SetInterval::interBox(const IntervalVector& box,NodeType valin,NodeType valout) {
+
+	// case root is leaf create issue for SetLeaf recursive function that need the father
+	// root is redefined as SetBisect pointing on two leaves with value of the root
+	if(root->is_leaf())
+		this->cutRoot();
+	cout<<"status of root: "<<root->status<<endl;
+	root->operator_ir(bounding_box,box,valin,valout,true,eps);
+	this->gather();
+	return *this;
+}
+
+SetInterval& SetInterval::unionBox(const IntervalVector& box,NodeType valin,NodeType valout) {
+	// case root is leaf create issue for SetLeaf recursive function that need the father
+	// root is redefined as SetBisect pointing on two leaves with value of the root
+	if(root->is_leaf())
+		this->cutRoot();
+	root->operator_ir(bounding_box,box,valin,valout,false,eps);
+	this->gather();
+	return *this;
+}
+
+
 SetInterval& SetInterval::operator|=(const SetInterval& set) {
-	root = root->union_(bounding_box, set.root, set.bounding_box, eps);
+	assert(bounding_box== set.bounding_box); // root must be the same as set interval are regular
+	// case root is leaf create issue for SetLeaf recursive function that need the father
+	// root is redefined as SetBisect pointing on two leaves with value of the root
+	if(root->is_leaf())
+		this->cutRoot();
+	root->oper(set.root,false);
+	this->gather();
 	return *this;
 }
 
+void SetInterval::gather() {
+	root->gather(false);
+	root->cutDeadBranch();
+}
+
+void SetInterval::checkFat() {
+	root->checkFat();
+}
+
 void SetInterval::save(const char* filename) {
 	std::stack<SetNode*> s;
 
@@ -175,6 +230,7 @@ void SetInterval::load(const char* filename) {
 	}
 
 	is.close();
+	root->setFathers();
 }
 
 void SetInterval::visit_leaves(SetNode::leaf_func func) const {
diff --git a/src/set/ibex_SetInterval.h b/src/set/ibex_SetInterval.h
index 120080c02..ad3982e6c 100644
--- a/src/set/ibex_SetInterval.h
+++ b/src/set/ibex_SetInterval.h
@@ -37,6 +37,12 @@ class SetInterval {
 	 */
 	SetInterval(const IntervalVector& bounding_box, double eps, bool inner=true);
 
+	/**
+	 * \brief Creates a copy of an existing setInterval
+	 *
+	 */
+	SetInterval(const SetInterval& set);
+
 
 	/**
 	 * \brief Loads a set from a data file.
@@ -50,6 +56,12 @@ class SetInterval {
 	 */
 	~SetInterval();
 
+
+	/*
+	 * \brief Replace the current root by a setbisect pointing on 2 leaves with the root status
+	 */	
+	void cutRoot();
+
 	/**
 	 * \brief i-Set Intersection
 	 *
@@ -60,6 +72,24 @@ class SetInterval {
 	 */
 	SetInterval& operator&=(const SetInterval& set);
 
+	/**
+	 * \brief i-Set Intersection
+	 *
+	 * Intersection of a setInterval with a box of interior value valin and exterior value valout. 
+	 * Computation is more efficient than create a setInterval contracted on the box, and then
+	 * compute intersection of the two setInterval. 
+	 */
+	SetInterval& interBox(const IntervalVector& box,NodeType valin,NodeType valout);
+
+	/**
+	 * \brief i-Set Intersection
+	 *
+	 * Union of a setInterval with a box of interior value valin and exterior value valout. 
+	 * Computation is more efficient than create a setInterval contracted on the box, and then
+	 * compute intersection of the two setInterval. 
+	 */
+	SetInterval& unionBox(const IntervalVector& box,NodeType valin,NodeType valout);
+
 	/**
 	 * \brief i-Set Union
 	 *
@@ -71,14 +101,14 @@ class SetInterval {
 	SetInterval& operator|=(const SetInterval& set);
 
 	/**
-	 * \brief i-Set synchronization
-	 *
-	 * In Jaulin's terminology, this operator is the "intersection of i-sets" (squared symbol)
-	 *
-	 * If [x] designates this i-set and [y] the i-set in argument, then this will be replace by
-	 *  { x, x\in[x] and x\in[y] }.
+	 * \brief gather leaves that have the same value
 	 */
-	void sync(Sep& sep);
+	void gather();
+
+	/**
+	 * \brief Check if fathers are initialized
+	 */
+	void checkFat();
 
 	/**
 	 * \brief True if this i-set is empty
@@ -87,6 +117,9 @@ class SetInterval {
 	 */
 	bool is_empty() const;
 
+	/**
+	 * \brief Contrat i-set w.r.t a separator sep
+	 */
 	void contract(Sep& sep);
 
 	/**
@@ -102,6 +135,10 @@ class SetInterval {
 	 */
 	double dist(const Vector& pt, bool inside) const;
 
+	
+
+	
+
 protected:
 
 	/**
diff --git a/src/set/ibex_SetLeaf.cpp b/src/set/ibex_SetLeaf.cpp
index f2386ee1a..a424c6b73 100644
--- a/src/set/ibex_SetLeaf.cpp
+++ b/src/set/ibex_SetLeaf.cpp
@@ -29,6 +29,11 @@ SetLeaf::SetLeaf(NodeType status) : SetNode(status) {
 		ibex_error("cannot set multiple status to SetLeaf");
 	}
 }
+
+SetLeaf::SetLeaf(NodeType status,SetNode* father): SetNode(status,father) {
+
+} 
+
 SetLeaf::~SetLeaf() {
 
 }
@@ -37,117 +42,282 @@ bool SetLeaf::is_leaf() const {
 	return true;
 }
 
-//SetNode* SetLeaf::sync(const IntervalVector& nodebox, const IntervalVector& x, NodeType xstatus, double eps, Mode mode) {
-//	switch (mode) {
-//	case SYNC: return sync(nodebox, x, xstatus, eps);
-//	case INTER: return inter(nodebox, x, xstatus, eps);
-//	case UNION: not_implemented("SetLeaf::sync with union"); return this; break;
-//	default : ibex_error("SetLeaf::sync: unknown mode"); return this;
-//	}
-//}
-
-SetNode* SetLeaf::sync(const IntervalVector& nodebox, const IntervalVector& x, NodeType xstatus, double eps) {
-	//cout << nodebox << " " << to_string(status)  << " sync " << x << " ";
-	assert(xstatus<=UNK);
-
-	if (xstatus==UNK || status==xstatus) {
-		//cout << "this\n";
-		return this;
-	} else if (nodebox.is_subset(x)) {
-		if (status!=UNK) throw NoSet();
-		status=xstatus;
-		//cout << "this\n";
-		return this;
-	} else {
-		if ((nodebox & x).is_flat()) {
-			//cout << "this\n";
-			return this;
+ SetNode * SetLeaf::copy() const {
+	return new SetLeaf(this->status); 
+}
+
+void SetLeaf::inter(NodeType x_status) {
+	status = inte(status,x_status); 
+}
+
+void SetLeaf::_union(NodeType x) {
+	status = uni(status,x);
+}
+
+void SetLeaf::oper(SetNode * node,bool op) {
+	if(node->is_leaf()) // if node is leaf, apply its value to leaves of this
+	{
+		if(op)
+			status = inte(status,node->status);
+		else
+			status = uni(status,node->status);
+	}
+	else{ // if this is leaf and node isn't
+		SetBisect * bfather = (SetBisect*) father;
+		if(bfather->left == this)
+		{
+			bfather->left = node->copy(); // copy the structure of node
+			bfather->left->father = father;
+			if(op && status != IN)
+				bfather->left->inter(status);  // apply status to leaves
+			else if(!op && status!= OUT)
+				bfather->left->_union(status);
+		}
+		else
+		{
+			bfather->right = node->copy();
+			bfather->right->father = father;
+			if(op && status != IN)
+				bfather->right->inter(status);
+			else if(!op && status!= OUT)
+				bfather->right->_union(status);
 		}
-		SetNode* new_node=diff(nodebox, x, status, xstatus, eps);
-		delete this; // warning: suicide, don't move it before previous line
-		//cout << "gives "; new_node->print(cout,nodebox,0);
-		return new_node;
+		delete this; // delete this as it is now a setBisect and no more a leaf
+
 	}
 }
 
-SetNode* SetLeaf::sync_rec(const IntervalVector& nodebox, Sep& sep, double eps) {
-
-	if (status<UNK || nodebox.max_diam()<=eps) {
-		return this;
-	} else {
-		int var=nodebox.extr_diam_index(false);
-		pair<IntervalVector,IntervalVector> p=nodebox.bisect(var);
-		double pt=p.first[var].ub();
-		assert(nodebox[var].interior_contains(pt));
-		SetBisect* bis = new SetBisect(var, pt, new SetLeaf(UNK), new SetLeaf(UNK));
-		delete this;
-		return bis->sync_rec(nodebox, sep, eps);
+void SetLeaf::operator_ir(const IntervalVector& box,const IntervalVector& subbox, NodeType val, bool op,double eps) {
+    if(box.is_subset(subbox))
+	{
+		if(op && val != IN)
+			status = inte(status,val);
+		else if(!op && val!=OUT)
+			status = uni(status,val);
+	}
+	else if(box.max_diam()>eps)
+	{
+		SetBisect * bfather = (SetBisect*) father;
+		int var = box.extr_diam_index(false); // indice of the maximal diameter
+		double pt = box[var].mid();
+		SetLeaf * lnode = new SetLeaf(status,NULL); // create future right node
+		SetLeaf * rnode = new SetLeaf(status,NULL); // create future left node
+        if(bfather->left == this) // this is father->left
+		{
+			bfather->left = new SetBisect(var,pt,lnode,rnode,father); // create SetBisect to replace leaf
+			lnode->father = bfather->left;  // set fathers
+			rnode->father = bfather->left;
+            bfather->left->operator_ir(box,subbox,val,op,eps);
+		}
+		else // this is father->right
+		{
+			bfather->right = new SetBisect(var,pt,lnode,rnode,father);
+			lnode->father = bfather->right;
+			rnode->father = bfather->right;
+            bfather->right->operator_ir(box,subbox,val,op,eps);
+		}
+		delete this; // delete former father->left or father->right
 	}
+	else // if leaf reach minimum precision, compute partial intersection or union
+	{
+		if(op)
+			status = inte_in(status,val);
+		else
+			status = uni_in(status,val);
+    }
 }
 
+// Uncomment this function and comment the one above to use fakeBranch method
+/*void SetLeaf::operator_ir(const IntervalVector& box,const IntervalVector& subbox, NodeType val, bool op,bool is_left,double eps) {
+	if(box.is_subset(subbox))
+	{
+		if(op)
+			status = inte(status,val);
+		else
+			status = uni(status,val);
+	}
+	else if(box.overlaps(subbox))
+	{
+		if(box.max_diam()>eps)
+		{
+			SetBisect * bfather = (SetBisect*) father;
+			int var = box.extr_diam_index(false); // indice of the maximal diameter
+			double pt = box[var].mid();
+			if(bfather->left == this)
+			{
+				bfather->left = new SetBisect(var,pt,NULL,NULL,father);
+				bfather->left->status = status; // create new set bissect with nodetype IN OUT or UNK, that it will transmit this value to its left and right
+				bfather->left->operator_ir(box,subbox,val,op,eps);
+			}
+			else
+			{
+				bfather->right = new SetBisect(var,pt,NULL,NULL,father);
+				bfather->right->status = status; // create new set bissect with nodetype IN OUT or UNK, that it will transmit this value to its left and right
+				bfather->right->operator_ir(box,subbox,val,op,eps);
+			}
+			delete this;
+		}
+		else
+			if(op)
+				status = inte_in(status,val);
+			else
+				status = uni_in(status,val);
+	}
+	else
+		return;
 
+}*/
 
-SetNode* SetLeaf::inter(const IntervalVector& nodebox, const IntervalVector& x, NodeType xstatus, double eps) {
-	//cout << nodebox << " " << to_string(status)  << " inter " << x << " ";
-	assert(xstatus<=UNK);
-
-	if (status<UNK || xstatus==UNK) {
-		//cout << "this\n";
-		return this;
-	} else if (nodebox.is_subset(x)) {
-		if (xstatus==IN && status==IN_TMP) status=IN; // if status==UNK, it remains UNK.
-		else if (xstatus==OUT ) status=OUT;
-		 //cout << "this\n";
-		return this;
-	} else if (status==UNK) {
-		 //cout << "this\n";
-		return this;
-	} else {
-		// status=(IN_TMP), xstatus=(IN | OUT).
-		SetNode* new_node=diff(nodebox, x, status, xstatus==OUT? OUT : IN, eps);
-		delete this; // warning: suicide, don't move it before previous line
-		//cout << "gives "; new_node->print(cout,nodebox,0);
-		return new_node;
+void SetLeaf::operator_ir(const IntervalVector& box,const IntervalVector& subbox, NodeType valin,NodeType valout, bool op,double eps) {
+    if(box.is_subset(subbox))
+	{
+		if(op && valin!=IN)
+			status = inte(status,valin);
+		else if(!op && valin!=OUT)
+			status = uni(status,valin);
 	}
+	else if(box.max_diam()>eps)
+	{
+		SetBisect * bfather = (SetBisect*) father;
+		int var = box.extr_diam_index(false); // indice of the maximal diameter
+		double pt = box[var].mid();
+		SetLeaf * lnode = new SetLeaf(status,NULL); // create future right node
+		SetLeaf * rnode = new SetLeaf(status,NULL); // create future left node
+        if(bfather->left == this) // this is father->left
+		{
+			bfather->left = new SetBisect(var,pt,lnode,rnode,father); // create SetBisect to replace leaf
+			lnode->father = bfather->left;  // set fathers
+			rnode->father = bfather->left;
+            bfather->left->operator_ir(box,subbox,valin,valout,op,eps);
+		}
+		else // this is father->right
+		{
+			bfather->right = new SetBisect(var,pt,lnode,rnode,father);
+			lnode->father = bfather->right;
+			rnode->father = bfather->right;
+            bfather->right->operator_ir(box,subbox,valin,valout,op,eps);
+		}
+		delete this; // delete former father->left or father->right
+	}
+	else // if leaf reach minimum precision, compute partial intersection or union
+	{
+		if(op)
+			status = inte_in(status,UNK);
+		else
+			status = uni_in(status,UNK);
+    }
 }
 
-SetNode* SetLeaf::inter_rec(const IntervalVector& nodebox, Sep& sep, double eps) {
-
-	if (status<UNK)
-		return this;
-	else if (nodebox.max_diam()<=eps) {
-		return this;  // status is IN_TMP and stays like this.
-	} else {
-		int var=nodebox.extr_diam_index(false);
-		pair<IntervalVector,IntervalVector> p=nodebox.bisect(var);
-		double pt=p.first[var].ub();
-		assert(nodebox[var].interior_contains(pt));
-		SetBisect* bis = new SetBisect(var, pt, new SetLeaf(status), new SetLeaf(status));
-		delete this;
-		return bis->inter_rec(nodebox, sep, eps);
-	}
+SetNode * SetLeaf::fakeLeaf(const IntervalVector& box,const IntervalVector& subbox,NodeType val,bool op, double eps) {
+	assert(!this->is_leaf());
+	return NULL;
 }
 
-SetNode* SetLeaf::union_(const IntervalVector& nodebox, const IntervalVector& x, NodeType xstatus, double eps) {
-	//cout << nodebox << " " << to_string(status)  << " union << x << " ";
-	assert(xstatus<=UNK);
-
-	if (status==IN || xstatus>IN) {
-		//cout << "this\n";
-		return this;
-	} else if (nodebox.is_subset(x)) {
-		status=IN;
-		 //cout << "this\n";
-		return this;
-	} else {
-		// status=(UNK | OUT), xstatus=(IN).
-		SetNode* new_node=diff(nodebox, x, status, IN, eps);
-		delete this; // warning: suicide, don't move it before previous line
-		//cout << "gives "; new_node->print(cout,nodebox,0);
-		return new_node;
+void SetLeaf::cleave(const IntervalVector& box, Sep& sep, const double eps) {
+
+	/*IntervalVector box1(box);
+	IntervalVector box2(box);
+	sep.separate(box1,box2);
+	if(box1.is_empty())
+		this->inter(OUT);
+	else if(box2.is_empty())
+		this->_union(IN);
+	else if(box.max_diam()>eps)
+	{
+		SetBisect * bfather = (SetBisect*) father;
+		int var = box.extr_diam_index(false); // indice of the maximal diameter
+		double pt = box[var].mid();
+		SetLeaf * lnode = new SetLeaf(status,NULL); // create future right node
+		SetLeaf * rnode = new SetLeaf(status,NULL); // create future left node
+        if(bfather->left == this) // this is father->left
+		{
+			bfather->left = new SetBisect(var,pt,lnode,rnode,father); // create SetBisect to replace leaf
+			lnode->father = bfather->left;  // set fathers
+			rnode->father = bfather->left;
+            bfather->left->cleave(box,sep,eps);
+		}
+		else // this is father->right
+		{
+			bfather->right = new SetBisect(var,pt,lnode,rnode,father);
+			lnode->father = bfather->right;
+			rnode->father = bfather->right;
+            bfather->right->cleave(box,sep,eps);
+		}
+		delete this; // delete former father->left or father->right
+	}
+	else
+		status = inte(status,UNK);*/
+
+	assert(status == UNK);
+	IntervalVector box1(box);
+	IntervalVector box2(box);
+	sep.separate(box1,box2);
+	if(box1.is_empty()) // all the box can be set to OUT
+	{
+		this->inter(OUT);
+		return;
+	}
+	else if(box2.is_empty())// all the box can be set to IN
+	{
+		this->_union(IN);
+		return;
 	}
+	else if(box.max_diam()>eps)
+	{
+		SetBisect * bfather = (SetBisect*) father;
+		IntervalVector* restout;
+		int nout = 0;
+		if(box1!=box)
+    		nout=box.diff(box1,restout);
+    	IntervalVector* restin;
+    	int nin = 0;
+    	if(box2!=box)
+    		nin=box.diff(box2,restin); 	
+    	
+        if(bfather->left == this) // this is father->left
+		{
+			for(int i = 0;i<nout;i++)
+    			bfather->left->operator_ir(box,restout[i],OUT,true,eps); // add in box
+    		for(int i = 0;i<nin;i++)
+    			bfather->left->operator_ir(box,restin[i],IN,false,eps); // add out box
+    		bfather->left->gatherTo(false,bfather->left);
+    		if(!bfather->left->is_leaf())
+    			bfather->left->cleave(box,sep,eps);
+		}
+		else // this is father->right
+		{
+			for(int i = 0;i<nout;i++)
+    			bfather->right->operator_ir(box,restout[i],OUT,true,eps); // add in box
+    		for(int i = 0;i<nin;i++)
+    			bfather->right->operator_ir(box,restin[i],IN,false,eps); // add out box
+    		bfather->right->gatherTo(false,bfather->right);
+    		if(!bfather->right->is_leaf())
+    			bfather->right->cleave(box,sep,eps);
+		}
+	}
+	else
+		status = inte(status,UNK);
+	
+
+}
+
+void SetLeaf::gather(bool go_up) {
+	return; // nothing to do
+}
+
+void SetLeaf::gatherTo(bool go_up,SetNode * branch) {
+	return; // nothing to do
+}
+
+void SetLeaf::cutDeadBranch() {
+	return; // nothing to do
+}
+void SetLeaf::checkFat() {
+	if (father == NULL)
+		cout<<"issue, invalid father"<<endl;
 }
 
+
 void SetLeaf::visit_leaves(leaf_func func, const IntervalVector& nodebox) const {
 	func(nodebox, status==IN? YES : (status==OUT? NO : MAYBE));
 }
@@ -157,12 +327,8 @@ void SetLeaf::print(ostream& os, const IntervalVector& nodebox, int shift) const
 	os  << nodebox << " " << to_string(status) << endl;
 }
 
-void SetLeaf::set_in_tmp() {
-	if (status==IN) status=IN_TMP;
-}
-
-void SetLeaf::unset_in_tmp() {
-	if (status==IN_TMP) status=UNK;
+void SetLeaf::setFathers() {
+	return; // nothing to do here
 }
-
 } // namespace ibex
+
diff --git a/src/set/ibex_SetLeaf.h b/src/set/ibex_SetLeaf.h
index b66a2c443..b5007dc62 100644
--- a/src/set/ibex_SetLeaf.h
+++ b/src/set/ibex_SetLeaf.h
@@ -26,30 +26,49 @@ class SetLeaf : public SetNode {
 	 */
 	SetLeaf(NodeType status);
 
+	SetLeaf(NodeType status,SetNode* father);
+
 	/**
 	 * \brief Delete this.
 	 */
 	virtual ~SetLeaf();
 
 	/** \see SetNode */
+		/** \see SetNode */
 	virtual bool is_leaf() const;
 
-	//virtual SetNode* sync(const IntervalVector& nodebox, const IntervalVector& x, NodeType x_status, double eps, Mode mode);
+	/** \see SetNode */
+	virtual SetNode * copy() const;
+
+	/** \see SetNode */
+	virtual void inter(NodeType x_status);
 
 	/** \see SetNode */
-	virtual SetNode* sync(const IntervalVector& nodebox, const IntervalVector& x, NodeType x_status, double eps);
+	virtual void _union(NodeType x);
 
 	/** \see SetNode */
-	virtual SetNode* sync_rec(const IntervalVector& nodebox, Sep& sep, double eps);
+	virtual void oper(SetNode * other,bool op);
 
 	/** \see SetNode */
-	virtual SetNode* inter(const IntervalVector& nodebox, const IntervalVector& x, NodeType x_status, double eps);
+	virtual void operator_ir(const IntervalVector& box,const IntervalVector& nodebox, NodeType val,bool op, double eps);
 
 	/** \see SetNode */
-	virtual SetNode* inter_rec(const IntervalVector& nodebox, Sep& sep, double eps);
+	virtual void operator_ir(const IntervalVector& box,const IntervalVector& subbox, NodeType valin,NodeType valout, bool op,double eps);
 
 	/** \see SetNode */
-	virtual SetNode* union_(const IntervalVector& nodebox, const IntervalVector& x, NodeType x_status, double eps);
+	virtual SetNode * fakeLeaf(const IntervalVector& box,const IntervalVector& subbox,NodeType val,bool op, double eps);
+
+	/** \see SetNode */
+	virtual void cleave(const IntervalVector& box, Sep& sep, double eps);
+
+	/** \see SetNode */
+	virtual void gather(bool go_up);
+
+	/** \see SetNode */
+	virtual void cutDeadBranch();
+
+	/** \see SetNode */
+	virtual void checkFat();
 
 	/** \see SetNode */
 	virtual void visit_leaves(leaf_func func, const IntervalVector& nodebox) const;
@@ -58,10 +77,8 @@ class SetLeaf : public SetNode {
 	virtual void print(std::ostream& os, const IntervalVector& nodebox, int shift) const;
 
 	/** \see SetNode */
-	virtual void set_in_tmp();
+	virtual void setFathers();
 
-	/** \see SetNode */
-	virtual void unset_in_tmp();
 
 private:
 	SetLeaf(const SetLeaf&); // forbidden
diff --git a/src/set/ibex_SetNode.cpp b/src/set/ibex_SetNode.cpp
index df1e7926a..35824f5f9 100644
--- a/src/set/ibex_SetNode.cpp
+++ b/src/set/ibex_SetNode.cpp
@@ -26,117 +26,6 @@ using namespace std;
 
 namespace ibex {
 
-SetNode* diff(const IntervalVector& x, const IntervalVector& y, NodeType x_status, NodeType y_status, double eps) {
-
-	if (y.is_empty()) {
-		return new SetLeaf(x_status);
-	}
-
-	const int nn=x.size();
-
-	typedef struct {
-		int var;
-		double pt;
-		bool y_left; // is the part that belongs to y on the left (=true) or right (=false) side?
-	} bisection;
-
-	bisection *tmp = new bisection[2*nn]; // in the worst case, there is 2n bisections
-	Interval c1, c2;
-	int b=0;
-
-	for (int var=0; var<nn; var++) {
-
-		x[var].diff(y[var],c1,c2);
-		//cout << "x[" << var << "]=" << x[var] << " y[" << var << "]=" << y[var] << " c1=" << c1 << " c2=" << c2 << endl;
-
-		if (c1.is_empty()) continue;
-
-		// we discard the x-part if it is not significant
-		// but this is safe only if the status of y is UNK.
-		if (y_status==UNK && c1.diam()<eps) continue;
-
-		// when the set difference is close to x itself
-		// we return x. But this is safe only if the
-		// status of x is UNK.
-		if (x[var]==c1 || (x_status==UNK && x[var].delta(c1)<eps)) {
-			// no (significant) difference at all
-			delete[] tmp;
-			return new SetLeaf(x_status);
-		}
-
-		tmp[b].var = var;
-		if (c1.lb()==x[var].lb()) {
-			tmp[b].pt = c1.ub();
-			assert(x[var].interior_contains(tmp[b].pt));
-			tmp[b].y_left = false;
-		} else {
-			tmp[b].pt = c1.lb();
-			assert(x[var].interior_contains(tmp[b].pt));
-			tmp[b].y_left = true;
-		}
-		b++;
-
-		if (c2.is_empty()) continue;
-		if (y_status==UNK && c2.diam()<eps)  continue;
-		if (x[var]==c2 || (x_status==UNK && x[var].delta(c2)<eps)) {
-			// no (significant) difference at all
-			delete[] tmp;
-			return new SetLeaf(x_status);
-		}
-
-		tmp[b].var = var;
-		tmp[b].pt = c2.lb();
-		assert(x[var].interior_contains(tmp[b].pt));
-		tmp[b].y_left = true;
-		b++;
-	}
-
-	SetNode* root;
-
-	if (b==0) {
-		root = new SetLeaf(y_status);
-	} else {
-		SetNode* bisect=new SetLeaf(y_status);
-
-		for (int i=b-1; i>=0; i--) {
-			if (tmp[i].y_left) {
-				assert(x[tmp[i].var].interior_contains(tmp[i].pt));
-				bisect=new SetBisect(tmp[i].var,tmp[i].pt, bisect, new SetLeaf(x_status));
-			} else {
-				assert(x[tmp[i].var].interior_contains(tmp[i].pt));
-				bisect=new SetBisect(tmp[i].var,tmp[i].pt, new SetLeaf(x_status), bisect);
-			}
-		}
-		root = bisect;
-	}
-	delete[] tmp;
-
-	return root;
-}
-
-// result of a contraction represented as a SetNode
-SetNode* contract_set(const IntervalVector& box, Sep& sep, double eps) {
-
-	IntervalVector box1(box);
-	IntervalVector box2(box);
-
-	sep.separate(box1,box2);
-
-	SetNode* root1 = diff(box, box1, IN, UNK, eps);
-	//cout << "set obtained with inner contraction:" << endl; root1->print(cout,box,0);
-	SetNode* root2 = diff(box, box2, OUT, UNK, eps);
-	//cout << "set obtained with outer contraction:" << endl;	root2->print(cout,box,0);
-
-	// TODO: pb: sync fails
-	// check "false" (skip_other_maybe is an old parameter not used anymore)
-	root1->inter(box,root2,box,eps); //,false);
-
-	delete root2;
-
-	//cout << "final set:" << endl; root1->print(cout,box,0);
-
-	return root1;
-}
 
 char to_string(const NodeType& status) {
 	switch(status) {
@@ -146,99 +35,22 @@ char to_string(const NodeType& status) {
 	}
 }
 
-SetNode::SetNode(NodeType status) : status(status) {
-
-}
-
-SetNode::~SetNode() {
-
-}
-
-SetNode* SetNode::sync(const IntervalVector& nodebox, Sep& sep, double eps) {
-	// perform contraction
-	//cout << "=== contract with ctc_in  =======" << endl;
-
-	// we skip other UNK-box if this node is not a leaf. This makes no difference
-	// if we are in SYNC mode but if we are in INTER mode, this prevents from
-	// this node to be "absorbed" by a temporary UNK box resulting from contraction.
-	SetNode* this2 = this->sync(nodebox, contract_set(nodebox, sep, eps), nodebox, eps, !is_leaf());
-
-	//cout << " sep gives: "; this2->print(cout,nodebox,0);
-	//cout << endl;
-
+SetNode::SetNode(NodeType status) : status(status), father(NULL) {
 
-	SetNode* this3 = this2->sync_rec(nodebox, sep, eps);
-
-	return this3;
 }
 
-SetNode* SetNode::sync(const IntervalVector& nodebox, const SetNode* other, const IntervalVector& otherbox, double eps, bool skip_other_maybe) {
-
-	if (nodebox.is_disjoint(otherbox))
-		return this;
-	else if (other->is_leaf()) {
-		//if (other->status<UNK) // || !skip_other_maybe)
-			return sync(nodebox, otherbox, other->status, eps);
-		//else
-		//	return this;
-	} else {
-
-		SetBisect* bisect_node = (SetBisect*) other;
-		SetNode* this2 = sync(nodebox, bisect_node->left, bisect_node->left_box(otherbox), eps, skip_other_maybe);
-		// warning: cannot use this anymore (use this2 instead)
-		return this2->sync(nodebox, bisect_node->right, bisect_node->right_box(otherbox), eps, skip_other_maybe);
-	}
-}
-
-SetNode* SetNode::inter(const IntervalVector& nodebox, Sep& sep, double eps) {
-	// we skip other UNK-box if this node is not a leaf. This makes no difference
-	// if we are in SYNC mode but if we are in INTER mode, this prevents from
-	// this node to be "absorbed" by a temporary UNK box resulting from contraction.
-	SetNode* this2 = this->inter(nodebox, contract_set(nodebox, sep, eps), nodebox, eps);
-	//cout << " sep gives: "; this2->print(cout,nodebox,0);
-
-	SetNode* this4 = this2->inter_rec(nodebox, sep, eps);
+SetNode::SetNode(NodeType status, SetNode *father): status(status), father(father) {
 
-	return this4;
 }
 
+SetNode::~SetNode() {
 
-SetNode* SetNode::inter(const IntervalVector& nodebox, const SetNode* other, const IntervalVector& otherbox, double eps) {
-
-	if (nodebox.is_disjoint(otherbox))
-		return this;
-	else if (other->is_leaf()) {
-		// we consider IN_TMP to be "UNK" until the end. Otherwise, there
-		// would be problems when other->status==UNK. Indeed, if we set the status
-		// of this node to UNK, we lose information ("other" may be an intermediate
-		// node in the process, although it is a SetLeaf, and some sub-nodes of other
-		// could be IN. But once a node is UNK it cannot be set to IN anymore...).
-		// If we impose this node to be a leaf, it does not fix the problem (this node may also be an
-		// intermediate node in the process: consider a plain box (a leaf) to be
-		// contracted).
-			return inter(nodebox, otherbox, other->status, eps);
-	} else {
-
-		SetBisect* bisect_node = (SetBisect*) other;
-		SetNode* this2 = inter(nodebox, bisect_node->left, bisect_node->left_box(otherbox), eps);
-		// warning: cannot use this anymore (use this2 instead)
-		return this2->inter(nodebox, bisect_node->right, bisect_node->right_box(otherbox), eps);
-	}
 }
 
-SetNode* SetNode::union_(const IntervalVector& nodebox, const SetNode* other, const IntervalVector& otherbox, double eps) {
+/*void SetNode::inter(const SetNode* other)
+{
 
-	if (nodebox.is_disjoint(otherbox))
-		return this;
-	else if (other->is_leaf()) {
-		return union_(nodebox, otherbox, other->status, eps);
-	} else {
+}*/
 
-		SetBisect* bisect_node = (SetBisect*) other;
-		SetNode* this2 = union_(nodebox, bisect_node->left, bisect_node->left_box(otherbox), eps);
-		// warning: cannot use this anymore (use this2 instead)
-		return this2->union_(nodebox, bisect_node->right, bisect_node->right_box(otherbox), eps);
-	}
-}
 
 } // namespace ibex
diff --git a/src/set/ibex_SetNode.h b/src/set/ibex_SetNode.h
index 09dd5fbfa..a6a289c83 100644
--- a/src/set/ibex_SetNode.h
+++ b/src/set/ibex_SetNode.h
@@ -26,6 +26,16 @@ namespace ibex {
 class NoSet {
 
 };
+/**
+ * \ingroup iset
+ * \brief Exception thrown by sync function.
+ *
+ * use to compute regular cutting of a interval
+ */
+class vcut{
+	int var;
+	double pt;
+};
 
 /**
  * \brief Set node.
@@ -43,6 +53,12 @@ class SetNode {
 	 */
 	SetNode(NodeType status);
 
+	/**
+	 * \brief Creates a node of given status and father
+	 */
+	 SetNode(NodeType status,SetNode *father);
+
+
 	/**
 	 * \brief Delete this.
 	 */
@@ -54,69 +70,64 @@ class SetNode {
 	virtual bool is_leaf() const=0;
 
 	/**
-	 * \brief Synchronization with an i-set represented implicitly by a Sep
+	 * \brief Return a copy of the node.
 	 */
-	SetNode* sync(const IntervalVector& nodebox, Sep& sep, double eps);
+	virtual SetNode * copy() const=0;
+
 
 	/**
-	 * \brief Synchronization with an explicit i-set "other"
-	 *
-	 * skip_other_maybe: don't consider UNK box of the other set. This is important
-	 *                   because the other set may be a temporary set produced by
-	 *                   the contract function (and the UNK box will be refined later)
+	 * \brief Intersection or Union between a box of NodeType val and a regular setInterval.
+	 * Box may not respect the regularity of the regular setInterval.
+	 * The part of the set that intersect nodebox only is modify.
 	 */
-	SetNode* sync(const IntervalVector& nodebox, const SetNode* other, const IntervalVector& otherbox, double eps, bool skip_other_maybe);
+	virtual void operator_ir(const IntervalVector& box,const IntervalVector& nodebox, NodeType val, bool op, double eps)=0;
 
 	/**
-	 * \brief Synchronization with an i-set reduced to a single box "x" of status "x_status".
-	 *
+	 * \brief Call by operator_ir, create a subtree from a leaf with "fake leves" that are setbisect with left and right set to null
+	 * should avoid to create leaf and eventually destruct it if it need to be cut i.e replace by a setbisect
+	 * method abandonned, oddly slower than creating leaves and replace them by setbisect.
 	 */
-	virtual SetNode* sync(const IntervalVector& nodebox, const IntervalVector& x, NodeType x_status, double eps)=0;
+	virtual SetNode * fakeLeaf(const IntervalVector& box,const IntervalVector& subbox,NodeType val,bool op, double eps)=0;
 
 	/**
-	 * \brief Synchronization with a Sep (recursive call)
-	 *
-	 * Once the current "nodebox" has been synchronized with the two sets obtained by applying the inner and
-	 * outer contraction, a recursive call is performed (where the Sep is applied on the sub-boxes).
-	 *
-	 * If this node is a leaf, it means it has to be bisected.
+	 * \brief Overloaded function, apply valout to the set outside of the box and valin to the inside, thus modify the whole set.
 	 */
-	virtual SetNode* sync_rec(const IntervalVector& nodebox, Sep& sep, double eps)=0;
+	virtual void operator_ir(const IntervalVector& box,const IntervalVector& subbox, NodeType valin,NodeType valout, bool op,double eps)=0;
 
 	/**
-	 * \brief Intersection with an i-set represented implicitly by a Sep
+	 * \brief Intersection between two setNode.
 	 */
-	SetNode* inter(const IntervalVector& nodebox, Sep& sep, double eps);
+	virtual void oper(SetNode* other,bool op)=0;
 
 	/**
-	 * \brief Intersection with an explicit i-set "other"
+	 * \brief Intersection between setnode and a status.
 	 */
-	SetNode* inter(const IntervalVector& nodebox, const SetNode* other, const IntervalVector& otherbox, double eps);
+	virtual void inter(NodeType x_status)=0;
 
 	/**
-	 * \brief Intersection with an i-set reduced to a single box "x" of status "x_status".
+	 * \brief Union between setnode and a status.
 	 */
-	virtual SetNode* inter(const IntervalVector& nodebox, const IntervalVector& x, NodeType x_status, double eps)=0;
+	virtual void _union(NodeType x)=0;
 
 	/**
-	 * \brief Intersection with a Sep (recursive call)
-	 *
-	 * Once the current "nodebox" has been intersected with the two sets obtained by applying the inner and
-	 * outer contraction, a recursive call is performed (where the Sep is applied on the sub-boxes).
-	 *
-	 * If this node is a leaf, it means it has to be bisected.
+	 * \brief Change value of branch if right and left got the same value IN, OUT or UNK.
 	 */
-	virtual SetNode* inter_rec(const IntervalVector& nodebox, Sep& sep, double eps)=0;
-
+	virtual void gather(bool go_up)=0;
 	/**
-	 * \brief Union with an explicit i-set "other"
+	 * \brief Delete the branch if left and right got the same value IN, OUT or UNK.
 	 */
-	SetNode* union_(const IntervalVector& nodebox, const SetNode* other, const IntervalVector& otherbox, double eps);
+	virtual void cutDeadBranch()=0;
 
 	/**
-	 * \brief Union with an i-set reduced to a single box "x" of status "x_status".
+	 * \brief Check if fathers are initialized.
 	 */
-	virtual SetNode* union_(const IntervalVector& nodebox, const IntervalVector& x, NodeType x_status, double eps)=0;
+	virtual void checkFat()=0;
+
+	/**
+	 * \brief Contrat i-set w.r.t separator sep, split leaves until boxin and boxout returned by sep are disjoints, 
+	 * and then call operator_ir to contract on them
+	 */	
+	virtual void cleave(const IntervalVector& box, Sep& sep,double eps)=0;
 
 	/**
 	 * \brief Visit the leaves of a tree with a callback "func"
@@ -129,36 +140,20 @@ class SetNode {
 	virtual void print(std::ostream& os, const IntervalVector& nodebox, int shift) const=0;
 
 	/**
-	 * Set all "IN" leaves to "IN_TMP".
-	 *
-	 * This function is useful for calculating set intersection only.
-	 * It sets all the inner nodes of the original set to IN_TMP before.
-	 * This allows to distinguish a node that is inside the original set (IN_TMP)
-	 * from a node that is proven to be inside the intersection (IN).
+	 * \brief Set fathers of the nodes
 	 */
-	virtual void set_in_tmp()=0;
+	virtual void setFathers()=0;
+
 
-	/**
-	 * Set all "IN_TMP" leaves to "UNK"
-	 *
-	 * This function is useful for calculating set intersection only.
-	 * Once the intersection is calculated, the nodes that was in the original
-	 * set but that could not be proven to be inside the intersection
-	 * have "UNK" status.
-	 */
-	virtual void unset_in_tmp()=0;
 
 	/**
 	 * \brief The status of the node
 	 */
 	NodeType status;
+	SetNode * father;
 };
 
 
-SetNode* diff(const IntervalVector& x, const IntervalVector& y, NodeType x_status, NodeType y_status, double eps);
-
-SetNode* contract_set(const IntervalVector& x, Sep& sep, double eps);
-
 } // namespace ibex
 
 #endif // __IBEX_SET_NODE_H__