diff --git a/src/basis_factorization/BasisFactorizationFactory.cpp b/src/basis_factorization/BasisFactorizationFactory.cpp
index 2352863c9..f87db4367 100644
--- a/src/basis_factorization/BasisFactorizationFactory.cpp
+++ b/src/basis_factorization/BasisFactorizationFactory.cpp
@@ -16,14 +16,14 @@
 #include "GlobalConfiguration.h"
 #include "LUFactorization.h"
 
-IBasisFactorization *BasisFactorizationFactory::createBasisFactorization( unsigned basisSize )
+IBasisFactorization *BasisFactorizationFactory::createBasisFactorization( unsigned basisSize, const IBasisFactorization::BasisColumnOracle &basisColumnOracle )
 {
     if ( GlobalConfiguration::BASIS_FACTORIZATION_TYPE ==
          GlobalConfiguration::LU_FACTORIZATION )
-        return new LUFactorization( basisSize );
+        return new LUFactorization( basisSize, basisColumnOracle );
     else if ( GlobalConfiguration::BASIS_FACTORIZATION_TYPE ==
               GlobalConfiguration::FORREST_TOMLIN_FACTORIZATION )
-        return new ForrestTomlinFactorization( basisSize );
+        return new ForrestTomlinFactorization( basisSize, basisColumnOracle );
 
     throw BasisFactorizationError( BasisFactorizationError::UNKNOWN_BASIS_FACTORIZATION_TYPE );
 }
diff --git a/src/basis_factorization/BasisFactorizationFactory.h b/src/basis_factorization/BasisFactorizationFactory.h
index 534a5be33..65c699b88 100644
--- a/src/basis_factorization/BasisFactorizationFactory.h
+++ b/src/basis_factorization/BasisFactorizationFactory.h
@@ -18,7 +18,7 @@
 class BasisFactorizationFactory
 {
 public:
-    static IBasisFactorization *createBasisFactorization( unsigned basisSize );
+    static IBasisFactorization *createBasisFactorization( unsigned basisSize, const IBasisFactorization::BasisColumnOracle &basisColumnOracle );
 };
 
 #endif // __BasisFactorizationFactory_h__
diff --git a/src/basis_factorization/ForrestTomlinFactorization.cpp b/src/basis_factorization/ForrestTomlinFactorization.cpp
index 3c0d98635..c0026168a 100644
--- a/src/basis_factorization/ForrestTomlinFactorization.cpp
+++ b/src/basis_factorization/ForrestTomlinFactorization.cpp
@@ -18,13 +18,13 @@
 #include <cstdlib>
 #include <cstring>
 
-ForrestTomlinFactorization::ForrestTomlinFactorization( unsigned m )
-    : _m( m )
+ForrestTomlinFactorization::ForrestTomlinFactorization( unsigned m, const BasisColumnOracle &basisColumnOracle )
+    : IBasisFactorization( basisColumnOracle )
+    , _m( m )
     , _B( NULL )
     , _Q( m )
     , _invQ( m )
     , _U( NULL )
-    , _invLP( NULL )
     , _explicitBasisAvailable( false )
     , _workMatrix( NULL )
     , _workVector( NULL )
@@ -43,11 +43,6 @@ ForrestTomlinFactorization::ForrestTomlinFactorization( unsigned m )
         throw BasisFactorizationError( BasisFactorizationError::ALLOCATION_FAILED,
                                        "ForrestTomlinFactorization::U" );
 
-    _invLP = new double[m * m];
-    if ( !_invLP )
-        throw BasisFactorizationError( BasisFactorizationError::ALLOCATION_FAILED,
-                                       "ForrestTomlinFactorization::invLPx" );
-
     for ( unsigned i = 0; i < _m; ++i )
     {
         _U[i] = new EtaMatrix( _m, i );
@@ -81,11 +76,6 @@ ForrestTomlinFactorization::ForrestTomlinFactorization( unsigned m )
     for ( unsigned row = 0; row < _m; ++row )
         _B[row * _m + row] = 1.0;
 
-    // Also, invLP = I
-    std::fill_n( _invLP, _m * _m, 0.0 );
-    for ( unsigned row = 0; row < _m; ++row )
-        _invLP[row * _m + row] = 1.0;
-
     // Us, Q and invQ are all initialized to the identity matrix anyway.
 }
 
@@ -112,12 +102,6 @@ ForrestTomlinFactorization::~ForrestTomlinFactorization()
 		_U = NULL;
 	}
 
-    if ( _invLP )
-    {
-        delete[] _invLP;
-        _invLP = NULL;
-    }
-
     if ( _workMatrix )
 	{
 		delete[] _workMatrix;
@@ -309,7 +293,7 @@ void ForrestTomlinFactorization::pushEtaMatrix( unsigned columnIndex, const doub
 
     // If the number of A matrices is too great, condense them.
     if ( _A.size() > GlobalConfiguration::REFACTORIZATION_THRESHOLD )
-        makeExplicitBasisAvailable();
+        refactorizeBasis();
 }
 
 void ForrestTomlinFactorization::forwardTransformation( const double *y, double *x ) const
@@ -533,8 +517,6 @@ void ForrestTomlinFactorization::storeFactorization( IBasisFactorization *other
     for ( const auto &lp : _LP )
         otherFTFactorization->_LP.append( lp->duplicate() );
 
-    memcpy( otherFTFactorization->_invLP, _invLP, sizeof(double) * _m * _m );
-
     List<AlmostIdentityMatrix *>::iterator aIt;
     for ( aIt = otherFTFactorization->_A.begin(); aIt != otherFTFactorization->_A.end(); ++aIt )
         delete *aIt;
@@ -572,8 +554,6 @@ void ForrestTomlinFactorization::restoreFactorization( const IBasisFactorization
     for ( const auto &lp : otherFTFactorization->_LP )
         _LP.append( lp->duplicate() );
 
-    memcpy( _invLP, otherFTFactorization->_invLP, sizeof(double) * _m * _m );
-
     List<AlmostIdentityMatrix *>::iterator aIt;
     for ( aIt = _A.begin(); aIt != _A.end(); ++aIt )
         delete *aIt;
@@ -691,56 +671,6 @@ void ForrestTomlinFactorization::initialLUFactorization()
         for ( unsigned j = 0; j <= i; ++j )
             u->_column[j] = _workMatrix[j * _m + i];
     }
-
-    // Compute inv(P1)inv(L1) ... inv(Ps)inv(LS) for later.
-
-    // Initialize to the identity matrix
-    std::fill_n( _invLP, _m * _m, 0 );
-    for ( unsigned i = 0; i < _m; ++i )
-        _invLP[i*_m + i] = 1.0;
-
-    // Multiply by inverted Ps and Ls
-    for ( auto lp = _LP.rbegin(); lp != _LP.rend(); ++lp )
-    {
-        if ( (*lp)->_pair )
-        {
-            // The inverse of a general permutation matrix is its transpose.
-            // However, since these are permutations of just two rows, the
-            // transpose is equal to the original - so no need to invert.
-            // Multiplying on the right permutes columns.
-            double temp;
-            unsigned first = (*lp)->_pair->first;
-            unsigned second = (*lp)->_pair->second;
-
-            for ( unsigned i = 0; i < _m; ++i )
-            {
-                temp = _invLP[i * _m + first];
-                _invLP[i * _m + first] = _invLP[i * _m + second];
-                _invLP[i * _m + second] = temp;
-            }
-		}
-        else
-        {
-            // The inverse of an eta matrix is an eta matrix with the special
-            // column negated and divided by the diagonal entry. The only
-            // exception is the diagonal entry itself, which is just the
-            // inverse of the original diagonal entry.
-            EtaMatrix *eta = (*lp)->_eta;
-            unsigned columnIndex = eta->_columnIndex;
-            double etaDiagonalEntry = 1 / eta->_column[columnIndex];
-
-            for ( unsigned row = 0; row < _m; ++row )
-            {
-                for ( unsigned col = columnIndex + 1; col < _m; ++col )
-                {
-                    _invLP[row * _m + columnIndex] -=
-                        ( eta->_column[col] * _invLP[row * _m + col] );
-                }
-
-                _invLP[row * _m + columnIndex] *= etaDiagonalEntry;
-            }
-        }
-    }
 }
 
 bool ForrestTomlinFactorization::explicitBasisAvailable() const
@@ -752,70 +682,8 @@ void ForrestTomlinFactorization::makeExplicitBasisAvailable()
 {
     if ( explicitBasisAvailable() )
         return;
-    /*
-      We know that the following equation holds:
-
-      Ak...A1 * LsPs...L1P1 * B = Q * Um...U1 * invQ
-
-      Or:
 
-      B = inv( Ak...A1 * LsPs...L1P1 ) * Q * Um...U1 * invQ
-        = inv(P1)inv(L1) ... inv(Ps)inv(LS) * inv(A1)...inv(Ak)
-          * Q * Um...U1 * invQ
-    */
-
-    // Start with the already-computed inv(P1)inv(L1) ... inv(Ps)inv(LS)
-    memcpy( _B, _invLP, sizeof(double) * _m * _m );
-
-    // Multiply by the inverse As. An inverse of an almost identity matrix
-    // is just that matrix with its special value negated (unless that
-    // matrix is diagonal.
-    for ( const auto &a : _A )
-    {
-        if ( a->_row == a->_column )
-        {
-            for ( unsigned j = 0; j < _m; ++j )
-                _B[j * _m + a->_column] *= ( 1 / a->_value );
-        }
-        else
-        {
-            for ( unsigned j = 0; j < _m; ++j )
-                _B[j * _m + a->_column] -= _B[j * _m + a->_row] * a->_value;
-        }
-    }
-
-    // Permute columns according to Q
-    for ( unsigned i = 0; i < _m; ++i )
-    {
-        unsigned columnToCopy = _invQ._ordering[i];
-        for ( unsigned j = 0; j < _m; ++j )
-            _workMatrix[j * _m + i] = _B[j * _m + columnToCopy];
-    }
-
-    // Multiply by Um...U1
-    for ( int i = _m - 1; i >= 0; --i )
-    {
-        for ( unsigned row = 0; row < _m; ++row )
-        {
-            for ( unsigned j = 0; j < _U[i]->_columnIndex; ++j )
-            {
-                _workMatrix[row * _m + _U[i]->_columnIndex] +=
-                    ( _U[i]->_column[j] * _workMatrix[row * _m + j] );
-            }
-        }
-    }
-
-    // Permute columns according to invQ
-    for ( unsigned i = 0; i < _m; ++i )
-    {
-        unsigned columnToCopy = _Q._ordering[i];
-        for ( unsigned j = 0; j < _m; ++j )
-            _B[j * _m + i] = _workMatrix[j * _m + columnToCopy];
-    }
-
-    clearFactorization();
-    initialLUFactorization();
-    _explicitBasisAvailable = true;
+    refactorizeBasis();
 }
 
 const double *ForrestTomlinFactorization::getBasis() const
@@ -988,12 +856,21 @@ void ForrestTomlinFactorization::dump() const
     printf( "\nDumping invQ:\n" );
     _invQ.dump();
 
-    printf( "*** Done dumping FT factorization ***\n\no" );
+    printf( "*** Done dumping FT factorization ***\n\n" );
 }
 
-const double *ForrestTomlinFactorization::getInvLP() const
+void ForrestTomlinFactorization::refactorizeBasis()
 {
-    return _invLP;
+    for ( unsigned column = 0; column < _m; ++column )
+    {
+        const double *basisColumn = _basisColumnOracle->getColumnOfBasis( column );
+        for ( unsigned row = 0; row < _m; ++row )
+            _B[row * _m + column] = basisColumn[row];
+    }
+
+    clearFactorization();
+    initialLUFactorization();
+    _explicitBasisAvailable = true;
 }
 
 //
diff --git a/src/basis_factorization/ForrestTomlinFactorization.h b/src/basis_factorization/ForrestTomlinFactorization.h
index 149e13175..f58fd3e87 100644
--- a/src/basis_factorization/ForrestTomlinFactorization.h
+++ b/src/basis_factorization/ForrestTomlinFactorization.h
@@ -33,7 +33,7 @@
 class ForrestTomlinFactorization : public IBasisFactorization
 {
 public:
-    ForrestTomlinFactorization( unsigned m );
+    ForrestTomlinFactorization( unsigned m, const BasisColumnOracle &basisColumnOracle );
     ~ForrestTomlinFactorization();
 
     /*
@@ -88,6 +88,11 @@ class ForrestTomlinFactorization : public IBasisFactorization
      */
     void invertBasis( double *result );
 
+    /*
+      Obtain the basis matrix from the oracle and compute a fresh factorization
+    */
+    void refactorizeBasis();
+
 public:
     /*
       For testing purposes only
@@ -97,7 +102,6 @@ class ForrestTomlinFactorization : public IBasisFactorization
     const EtaMatrix **getU() const;
     const List<LPElement *> *getLP() const;
     const List<AlmostIdentityMatrix *> *getA() const;
-    const double *getInvLP() const;
 
     void pushA( const AlmostIdentityMatrix &matrix );
     void setQ( const PermutationMatrix &Q );
@@ -124,13 +128,6 @@ class ForrestTomlinFactorization : public IBasisFactorization
     // A1 is the first element of the list
     List<AlmostIdentityMatrix *>_A;
 
-    /*
-      A helper matrix to facilitate computing the explicit basis.
-      This matrix stores:
-      inv(P1)inv(L1) ... inv(Ps)inv(LS)
-    */
-    double *_invLP;
-
     /*
       Indicates whether the explicit basis matrix is available.
     */
diff --git a/src/basis_factorization/IBasisFactorization.h b/src/basis_factorization/IBasisFactorization.h
index fc0f90ca8..d6d7b1d66 100644
--- a/src/basis_factorization/IBasisFactorization.h
+++ b/src/basis_factorization/IBasisFactorization.h
@@ -19,8 +19,19 @@ class IBasisFactorization
       This is the interfact class for a basis factorization.
     */
 public:
-    IBasisFactorization()
+    /*
+      A callback for obtaining columns of the basis matrix
+    */
+    class BasisColumnOracle
+    {
+    public:
+        virtual ~BasisColumnOracle() {}
+        virtual const double *getColumnOfBasis( unsigned column ) const = 0;
+    };
+
+    IBasisFactorization( const BasisColumnOracle &basisColumnOracle )
         : _factorizationEnabled( true )
+        , _basisColumnOracle( &basisColumnOracle )
     {
     }
 
@@ -96,6 +107,9 @@ class IBasisFactorization
       disabled.
     */
     bool _factorizationEnabled;
+
+protected:
+    const BasisColumnOracle *_basisColumnOracle;
 };
 
 #endif // __IBasisFactorization_h__
diff --git a/src/basis_factorization/LUFactorization.cpp b/src/basis_factorization/LUFactorization.cpp
index 507dd448c..108d7fafe 100644
--- a/src/basis_factorization/LUFactorization.cpp
+++ b/src/basis_factorization/LUFactorization.cpp
@@ -21,8 +21,9 @@
 #include "MStringf.h"
 #include "MalformedBasisException.h"
 
-LUFactorization::LUFactorization( unsigned m )
-    : _B0( NULL )
+LUFactorization::LUFactorization( unsigned m, const BasisColumnOracle &basisColumnOracle )
+    : IBasisFactorization( basisColumnOracle )
+    , _B0( NULL )
 	, _m( m )
     , _U( NULL )
     , _tempY( NULL )
@@ -121,9 +122,8 @@ void LUFactorization::pushEtaMatrix( unsigned columnIndex, const double *column
 
 	if ( ( _etas.size() > GlobalConfiguration::REFACTORIZATION_THRESHOLD ) && factorizationEnabled() )
 	{
-        log( "Number of etas exceeds threshold. Condensing and refactoring\n" );
-		condenseEtas();
-		factorizeMatrix( _B0 );
+        log( "Number of etas exceeds threshold. Refactoring basis\n" );
+        refactorizeBasis();
 	}
 }
 
@@ -158,41 +158,10 @@ void LUFactorization::LMultiplyLeft( const EtaMatrix *L, double *X ) const
 void LUFactorization::setBasis( const double *B )
 {
 	memcpy( _B0, B, sizeof(double) * _m * _m );
+    clearFactorization();
 	factorizeMatrix( _B0 );
 }
 
-void LUFactorization::condenseEtas()
-{
-    // Multiplication by an eta matrix on the right only changes one
-    // column of B0. The new column is a linear combination of the
-    // existing columns of B0, according to the eta column. We perform
-    // the computation in place.
-    for ( const auto &eta : _etas )
-    {
-		unsigned col = eta->_columnIndex;
-
-        // Iterate over the rows of B0
-		for ( unsigned i = 0; i < _m; ++i )
-        {
-            // Compute the linear combinations of the columns
-			double sum = 0.0;
-			for ( unsigned j = 0; j < _m; ++j )
-				sum += _B0[i * _m + j] * eta->_column[j];
-
-            if ( FloatUtils::isZero( sum ) )
-                sum = 0.0;
-
-            _B0[i * _m + col] = sum;
-		}
-
-		delete eta;
-	}
-
-    _etas.clear();
-    clearLPU();
-    factorizeMatrix( _B0 );
-}
-
 void LUFactorization::forwardTransformation( const double *y, double *x ) const
 {
     // If there's no LP factorization, it is implied that B0 = I.
@@ -356,7 +325,7 @@ void LUFactorization::rowSwap( unsigned rowOne, unsigned rowTwo, double *matrix
     }
 }
 
-void LUFactorization::clearLPU()
+void LUFactorization::clearFactorization()
 {
     List<LPElement *>::iterator element;
     for ( element = _LP.begin(); element != _LP.end(); ++element )
@@ -364,12 +333,13 @@ void LUFactorization::clearLPU()
     _LP.clear();
 
 	std::fill_n( _U, _m*_m, 0 );
+
+    _etas.clear();
 }
 
 void LUFactorization::factorizeMatrix( double *matrix )
 {
-    // Clear any previous factorization, initialize U
-	clearLPU();
+    // Initialize U
 	memcpy( _U, matrix, sizeof(double) * _m * _m );
 
 	for ( unsigned i = 0; i < _m; ++i )
@@ -447,12 +417,15 @@ void LUFactorization::storeFactorization( IBasisFactorization *other )
     ASSERT( _m == otherLUFactorization->_m );
     ASSERT( otherLUFactorization->_etas.size() == 0 );
 
-    // In order to reduce space requirements, condense the etas before storing a factorization
-    condenseEtas();
-    factorizeMatrix( _B0 );
+    // In order to reduce space requirements (for the etas), compute a fresh refactorization
+    if ( !_etas.empty() )
+        refactorizeBasis();
 
-    // Now we simply store _B0
-    otherLUFactorization->setBasis( _B0 );
+    // Store the new basis and factorization
+    memcpy( otherLUFactorization->_B0, _B0, sizeof(double) * _m * _m );
+    memcpy( otherLUFactorization->_U, _U, sizeof(double) * _m * _m );
+    for ( const auto &it : _LP )
+        otherLUFactorization->_LP.append( it->duplicate() );
 }
 
 void LUFactorization::restoreFactorization( const IBasisFactorization *other )
@@ -463,14 +436,13 @@ void LUFactorization::restoreFactorization( const IBasisFactorization *other )
     ASSERT( otherLUFactorization->_etas.size() == 0 );
 
     // Clear any existing data
-    for ( const auto &it : _etas )
-        delete it;
-
-    _etas.clear();
-	clearLPU();
+    clearFactorization();
 
-    // Store the new B0 and LU-factorize it
-    setBasis( otherLUFactorization->_B0 );
+    // Store the new basis and factorization
+    memcpy( _B0, otherLUFactorization->_B0, sizeof(double) * _m * _m );
+    memcpy( _U, otherLUFactorization->_U, sizeof(double) * _m * _m );
+    for ( const auto &it : otherLUFactorization->_LP )
+        _LP.append( it->duplicate() );
 }
 
 void LUFactorization::invertBasis( double *result )
@@ -553,7 +525,56 @@ bool LUFactorization::explicitBasisAvailable() const
 
 void LUFactorization::makeExplicitBasisAvailable()
 {
-    condenseEtas();
+    refactorizeBasis();
+}
+
+void LUFactorization::dump() const
+{
+    printf( "*** Dumping LU factorization ***\n\n" );
+
+    printf( "\nDumping LPs:\n" );
+    unsigned count = 0;
+    for ( const auto &lp : _LP )
+    {
+        printf( "LP[%i]:\n", _LP.size() - 1 - count );
+        ++count;
+        lp->dump();
+    }
+    printf( "\n\n" );
+
+    printf( "Dumping Us:\n" );
+    for ( unsigned col = 0; col < _m; ++col )
+    {
+        printf( "U[%u]:\n", col );
+        for ( unsigned i = 0; i < _m; ++i )
+        {
+            printf( "\t%lf\n", _U[i * _m + col] );
+        }
+        printf( "\n" );
+    }
+
+    printf( "Dumping etas:\n" );
+    for ( const auto &eta : _etas )
+    {
+        eta->dump();
+        printf( "\n" );
+    }
+
+    printf( "*** Done dumping LU factorization ***\n\n" );
+}
+
+void LUFactorization::refactorizeBasis()
+{
+    clearFactorization();
+
+    for ( unsigned column = 0; column < _m; ++column )
+    {
+        const double *basisColumn = _basisColumnOracle->getColumnOfBasis( column );
+        for ( unsigned row = 0; row < _m; ++row )
+            _B0[row * _m + column] = basisColumn[row];
+    }
+
+    factorizeMatrix( _B0 );
 }
 
 //
diff --git a/src/basis_factorization/LUFactorization.h b/src/basis_factorization/LUFactorization.h
index cb2d05a0a..53f8e9df0 100644
--- a/src/basis_factorization/LUFactorization.h
+++ b/src/basis_factorization/LUFactorization.h
@@ -25,7 +25,7 @@ class LPElement;
 class LUFactorization : public IBasisFactorization
 {
 public:
-    LUFactorization( unsigned m );
+    LUFactorization( unsigned m, const BasisColumnOracle &basisColumnOracle );
     ~LUFactorization();
 
     /*
@@ -107,12 +107,6 @@ class LUFactorization : public IBasisFactorization
     */
     void rowSwap( unsigned rowOne, unsigned rowTwo, double *matrix );
 
-    /*
-      Compute B0 * E1 ... *En for all stored eta matrices, and place
-      the result in B0.
-    */
-	void condenseEtas();
-
     /*
       Return true iff the basis matrix B0 is explicitly available.
     */
@@ -122,6 +116,7 @@ class LUFactorization : public IBasisFactorization
       Make the basis explicitly available
     */
     void makeExplicitBasisAvailable();
+    void refactorizeBasis();
 
     /*
       Get the explicit basis matrix
@@ -147,6 +142,11 @@ class LUFactorization : public IBasisFactorization
 	const List<LPElement *> getLP() const;
 	const List<EtaMatrix *> getEtas() const;
 
+    /*
+      Debug
+    */
+    void dump() const;
+
 private:
     /*
       The Basis matrix.
@@ -180,7 +180,7 @@ class LUFactorization : public IBasisFactorization
     /*
       Clear a previous factorization.
     */
-	void clearLPU();
+	void clearFactorization();
 
     /*
       Helper functions for backward- and forward-transformations.
diff --git a/src/basis_factorization/tests/MockColumnOracle.h b/src/basis_factorization/tests/MockColumnOracle.h
new file mode 100644
index 000000000..7fec144fa
--- /dev/null
+++ b/src/basis_factorization/tests/MockColumnOracle.h
@@ -0,0 +1,66 @@
+/*********************                                                        */
+/*! \file MockColumnOracle.h
+** \verbatim
+** Top contributors (to current version):
+**   Guy Katz
+** This file is part of the Marabou project.
+** Copyright (c) 2016-2017 by the authors listed in the file AUTHORS
+** in the top-level source directory) and their institutional affiliations.
+** All rights reserved. See the file COPYING in the top-level source
+** directory for licensing information.\endverbatim
+**/
+
+#ifndef __MockColumnOracle_h__
+#define __MockColumnOracle_h__
+
+#include "IBasisFactorization.h"
+
+class MockColumnOracle : public IBasisFactorization::BasisColumnOracle
+{
+public:
+    MockColumnOracle()
+        : _basis( NULL )
+    {
+    }
+
+    ~MockColumnOracle()
+    {
+        if ( _basis )
+        {
+            delete[] _basis;
+            _basis = NULL;
+        }
+    }
+
+    void storeBasis( unsigned m, const double *basis )
+    {
+        _m = m;
+        _basis = new double[_m * _m];
+
+        for ( unsigned row = 0; row < _m; ++row )
+        {
+            for ( unsigned column = 0; column < _m; ++column )
+            {
+                _basis[column * _m + row] = basis[row * _m + column];
+            }
+        }
+    }
+
+    double *_basis;
+    unsigned _m;
+
+    const double *getColumnOfBasis( unsigned column ) const
+    {
+        return _basis + ( _m * column );
+    }
+};
+
+#endif // __MockColumnOracle_h__
+
+//
+// Local Variables:
+// compile-command: "make -C ../../.. "
+// tags-file-name: "../../../TAGS"
+// c-basic-offset: 4
+// End:
+//
diff --git a/src/basis_factorization/tests/Test_CompareFactorizations.h b/src/basis_factorization/tests/Test_CompareFactorizations.h
index 4ada62db3..34fa2dd0e 100644
--- a/src/basis_factorization/tests/Test_CompareFactorizations.h
+++ b/src/basis_factorization/tests/Test_CompareFactorizations.h
@@ -16,6 +16,7 @@
 #include "FloatUtils.h"
 #include "ForrestTomlinFactorization.h"
 #include "LUFactorization.h"
+#include "MockColumnOracle.h"
 #include "MockErrno.h"
 
 class MockForCompareFactorizations
@@ -27,14 +28,17 @@ class CompareFactorizationsTestSuite : public CxxTest::TestSuite
 {
 public:
     MockForCompareFactorizations *mock;
+    MockColumnOracle *oracle;
 
     void setUp()
     {
         TS_ASSERT( mock = new MockForCompareFactorizations );
+        TS_ASSERT( oracle = new MockColumnOracle );
     }
 
     void tearDown()
     {
+        TS_ASSERT_THROWS_NOTHING( delete oracle );
         TS_ASSERT_THROWS_NOTHING( delete mock );
     }
 
@@ -43,8 +47,8 @@ class CompareFactorizationsTestSuite : public CxxTest::TestSuite
         ForrestTomlinFactorization *ft;
         LUFactorization *lu;
 
-        TS_ASSERT( ft = new ForrestTomlinFactorization( 4 ) );
-        TS_ASSERT( lu = new LUFactorization( 4 ) );
+        TS_ASSERT( ft = new ForrestTomlinFactorization( 4, *oracle ) );
+        TS_ASSERT( lu = new LUFactorization( 4, *oracle ) );
 
         double B[] =
             {
@@ -87,7 +91,18 @@ class CompareFactorizationsTestSuite : public CxxTest::TestSuite
         for ( unsigned i = 0; i < 4; ++i )
             TS_ASSERT( FloatUtils::areEqual( x1[i], x2[i] ) );
 
+
+        double basisAtThisPoint[] = {
+            2, -20, 64, -4,
+            0, 2, 96, 0,
+            -3, 24, -45, 1,
+            0, 6, 14, 2,
+        };
+
+        oracle->storeBasis( 4, basisAtThisPoint );
+
         ft->makeExplicitBasisAvailable();
+        lu->makeExplicitBasisAvailable();
 
         TS_ASSERT_THROWS_NOTHING( ft->forwardTransformation( y, x1 ) );
         TS_ASSERT_THROWS_NOTHING( lu->forwardTransformation( y, x2 ) );
diff --git a/src/basis_factorization/tests/Test_ForrestTomlinFactorization.h b/src/basis_factorization/tests/Test_ForrestTomlinFactorization.h
index 5391e717b..d68242dd2 100644
--- a/src/basis_factorization/tests/Test_ForrestTomlinFactorization.h
+++ b/src/basis_factorization/tests/Test_ForrestTomlinFactorization.h
@@ -16,9 +16,10 @@
 #include "BasisFactorizationError.h"
 #include "EtaMatrix.h"
 #include "FloatUtils.h"
-#include "GlobalConfiguration.h"
 #include "ForrestTomlinFactorization.h"
+#include "GlobalConfiguration.h"
 #include "List.h"
+#include "MockColumnOracle.h"
 #include "MockErrno.h"
 
 class MockForForrestTomlinFactorization
@@ -30,6 +31,7 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite
 {
 public:
     MockForForrestTomlinFactorization *mock;
+    MockColumnOracle *oracle;
 
     bool isIdentityPermutation( const PermutationMatrix *matrix )
     {
@@ -43,10 +45,12 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite
     void setUp()
     {
         TS_ASSERT( mock = new MockForForrestTomlinFactorization );
+        TS_ASSERT( oracle = new MockColumnOracle );
     }
 
     void tearDown()
     {
+        TS_ASSERT_THROWS_NOTHING( delete oracle );
         TS_ASSERT_THROWS_NOTHING( delete mock );
     }
 
@@ -54,7 +58,7 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite
     {
         ForrestTomlinFactorization *ft;
 
-        TS_ASSERT( ft = new ForrestTomlinFactorization( 3 ) );
+        TS_ASSERT( ft = new ForrestTomlinFactorization( 3, *oracle ) );
 
         TS_ASSERT( ft->factorizationEnabled() );
 
@@ -73,7 +77,7 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite
     {
         ForrestTomlinFactorization *ft;
 
-        TS_ASSERT( ft = new ForrestTomlinFactorization( 4 ) );
+        TS_ASSERT( ft = new ForrestTomlinFactorization( 4, *oracle ) );
 
         double basisMatrix[16] = {
             1,   3, -2,  4,
@@ -188,48 +192,6 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite
         ++lIt;
         TS_ASSERT_EQUALS( *((*lIt)->_eta), expectedL1 );
 
-        /*
-          L4 = | 1 0 0    0 |
-               | 0 1 0    0 |
-               | 0 0 1    0 |
-               | 0 0 0 -0.5 |
-
-          L3 = | 1 0  0 0 |
-               | 0 1  0 0 |
-               | 0 0 -1 0 |
-               | 0 0  1 1 |
-
-          L2 = | 1 0   0 0 |
-               | 0 0.5 0 0 |
-               | 0 0   1 0 |
-               | 0 0   0 1 |
-
-          L1 = |  1 0 0 0 |
-               | -1 1 0 0 |
-               | -1 0 1 0 |
-               |  1 0 0 1 |
-
-         InvLP = inv( L4L3L2L1 )
-
-               = inv( | 1    0      0    0 | )
-                      | -0.5 0.5    0    0 |
-                      | 1    0     -1    0 |
-                      | 0    0   -0.5 -0.5 |
-
-               = | 1  0  0  0 |
-                 | 1  2  0  0 |
-                 | 1  0 -1  0 |
-                 | -1 0  1 -2 |
-        */
-        double expectedInvLP[] =
-            { 1, 0, 0, 0,
-              1, 2, 0, 0,
-              1, 0, -1, 0,
-              -1, 0, 1, - 2};
-
-        for ( unsigned i = 0; i < 16; ++i )
-            TS_ASSERT( FloatUtils::areEqual( expectedInvLP[i], ft->getInvLP()[i] ) );
-
         TS_ASSERT_THROWS_NOTHING( delete ft );
     }
 
@@ -237,7 +199,7 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite
     {
         ForrestTomlinFactorization *ft;
 
-        TS_ASSERT( ft = new ForrestTomlinFactorization( 4 ) );
+        TS_ASSERT( ft = new ForrestTomlinFactorization( 4, *oracle ) );
 
         double basisMatrix[16] = {
             1,   4,  -2,  4,
@@ -373,7 +335,7 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite
     {
         ForrestTomlinFactorization *ft;
 
-        TS_ASSERT( ft = new ForrestTomlinFactorization( 4 ) );
+        TS_ASSERT( ft = new ForrestTomlinFactorization( 4, *oracle ) );
 
         double basisMatrix[16] = {
             1,   3, -2,  4,
@@ -501,7 +463,7 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite
     {
         ForrestTomlinFactorization *ft;
 
-        TS_ASSERT( ft = new ForrestTomlinFactorization( 4 ) );
+        TS_ASSERT( ft = new ForrestTomlinFactorization( 4, *oracle ) );
 
         double basisMatrix[16] = {
             1,   3, -2,  4,
@@ -640,7 +602,7 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite
     {
         ForrestTomlinFactorization *ft;
 
-        TS_ASSERT( ft = new ForrestTomlinFactorization( 4 ) );
+        TS_ASSERT( ft = new ForrestTomlinFactorization( 4, *oracle ) );
 
         // B = | 1   3 -2  4 |
         //     | 1   5 -1  5 |
@@ -854,7 +816,7 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite
     {
         ForrestTomlinFactorization *ft;
 
-        TS_ASSERT( ft = new ForrestTomlinFactorization( 4 ) );
+        TS_ASSERT( ft = new ForrestTomlinFactorization( 4, *oracle ) );
 
         double basisMatrix[16] = {
             1,   3, -2,  4,
@@ -868,8 +830,8 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite
         double a1[] = { -4, 2, 0, 3 };
         ft->pushEtaMatrix( 1, a1 );
 
-        ForrestTomlinFactorization *ft2 = new ForrestTomlinFactorization( 4 );
-        ForrestTomlinFactorization *ft3 = new ForrestTomlinFactorization( 4 );
+        ForrestTomlinFactorization *ft2 = new ForrestTomlinFactorization( 4, *oracle );
+        ForrestTomlinFactorization *ft3 = new ForrestTomlinFactorization( 4, *oracle );
 
         ft->storeFactorization( ft2 );
         ft3->restoreFactorization( ft2 );
@@ -894,7 +856,7 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite
     {
         ForrestTomlinFactorization *ft;
 
-        TS_ASSERT( ft = new ForrestTomlinFactorization( 4 ) );
+        TS_ASSERT( ft = new ForrestTomlinFactorization( 4, *oracle ) );
 
         // B = | 1   3 -2  4 |
         //     | 1   5 -1  5 |
@@ -934,6 +896,8 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite
             -1, -26,  3, -8,
         };
 
+        oracle->storeBasis( 4, expectedB );
+
         TS_ASSERT( !ft->explicitBasisAvailable() );
         TS_ASSERT_THROWS_NOTHING( ft->makeExplicitBasisAvailable() );
         TS_ASSERT( ft->explicitBasisAvailable() );
@@ -951,7 +915,7 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite
     {
         ForrestTomlinFactorization *ft;
 
-        TS_ASSERT( ft = new ForrestTomlinFactorization( 4 ) );
+        TS_ASSERT( ft = new ForrestTomlinFactorization( 4, *oracle ) );
 
         // B = | 1   3 -2  4 |
         //     | 1   5 -1  5 |
@@ -995,6 +959,15 @@ class ForrestTomlinFactorizationTestSuite : public CxxTest::TestSuite
         //          | 1   20 -3  6 |
         //          | -1 -26  3 -8 |
 
+        double expectedB[16] = {
+            1,   14, -2,  4,
+            1,   21, -1,  5,
+            1,   20, -3,  6,
+            -1, -26,  3, -8,
+        };
+
+        oracle->storeBasis( 4, expectedB );
+
         // invB = |  4   -1/2 -13/4  -3/4 |
         //        | -1/2  1/4   5/8   3/8 |
         //        |  1      0    -2    -1 |
diff --git a/src/basis_factorization/tests/Test_LUFactorization.h b/src/basis_factorization/tests/Test_LUFactorization.h
index 182cdd910..6a2cab76f 100644
--- a/src/basis_factorization/tests/Test_LUFactorization.h
+++ b/src/basis_factorization/tests/Test_LUFactorization.h
@@ -19,6 +19,7 @@
 #include "GlobalConfiguration.h"
 #include "LUFactorization.h"
 #include "List.h"
+#include "MockColumnOracle.h"
 #include "MockErrno.h"
 
 void matrixMultiply( unsigned dimension, const double *left, const double *right, double *result )
@@ -45,14 +46,17 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite
 {
 public:
     MockForLUFactorization *mock;
+    MockColumnOracle *oracle;
 
     void setUp()
     {
         TS_ASSERT( mock = new MockForLUFactorization );
+        TS_ASSERT( oracle = new MockColumnOracle );
     }
 
     void tearDown()
     {
+        TS_ASSERT_THROWS_NOTHING( delete oracle );
         TS_ASSERT_THROWS_NOTHING( delete mock );
     }
 
@@ -60,7 +64,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite
     {
         LUFactorization *basis;
 
-        TS_ASSERT( basis = new LUFactorization( 3 ) );
+        TS_ASSERT( basis = new LUFactorization( 3, *oracle ) );
 
         TS_ASSERT( basis->factorizationEnabled() );
 
@@ -77,7 +81,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite
 
     void test_forward_transformation()
     {
-        LUFactorization basis( 3 );
+        LUFactorization basis( 3, *oracle );
 
         // If no eta matrices are provided, d = a
         double a1[] = { 1, 1, 3 };
@@ -127,7 +131,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite
 	void test_forward_transformation_with_B0()
 	{
         // Same etas as test_backward_transformation()
-		LUFactorization basis( 3 );
+		LUFactorization basis( 3, *oracle );
 		double e1[] = {1., 1., 3.};
 		basis.pushEtaMatrix( 1, e1 );
 		double e2[] = {2., 1., 1.};
@@ -142,7 +146,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite
 
 		double a[] = { 2., -1., 4. };
 		double d[] = { 0., 0., 0. };
-        double expected[] = { 42, 116, -131 };
+        double expected[] = { -20, 27, -8 };
 
 		basis.forwardTransformation( a, d );
 
@@ -152,7 +156,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite
 
 	void test_backward_transformation()
     {
-        LUFactorization basis( 3 );
+        LUFactorization basis( 3, *oracle );
 
         // If no eta matrices are provided, x = y
         double y1[] = { 1, 2, 3 };
@@ -222,7 +226,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite
 
     void test_backward_transformation_2()
     {
-        LUFactorization basis( 3 );
+        LUFactorization basis( 3, *oracle );
 
         // E1 = | -1     |
         //      |  0 1   |
@@ -247,7 +251,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite
 	void test_backward_transformation_with_B0()
 	{
         // Same etas as test_backward_transformation()
-		LUFactorization basis( 3 );
+		LUFactorization basis( 3, *oracle );
 		double e1[] = {1., 1., 3.};
 		basis.pushEtaMatrix( 1, e1 );
 		double e2[] = {2., 1., 1.};
@@ -260,12 +264,13 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite
 
 		double y[] = {19., 12., 17.};
 		double x[] = {0., 0., 0.};
-		double expected[] = {-6, 9, -4};
-		//     	| 1 2 4	|  	| 1 1   |   | 2     |   | 1   0.5 |
-        //  x *	| 4	5 7 | * |   1   | * | 1	1	| *	|	1 0.5 | = | 19 12 17 |
-        //     	| 7 8 9	|	|   3 1 |   | 1   1 |   |     0.5 |
+		double expected[] = { -104.0/3, 140.0/3, -19 };
+
+		//     	| 1 2 4	|
+        //  x *	| 4	5 7 | = | 19 12 17 |
+        //     	| 7 8 9	|
         //
-        // --> x = [ -6 9 -4 ]
+        // --> x = [ -104/3, 140/3, -19 ]
 		basis.backwardTransformation( y, x );
 
         for ( unsigned i = 0; i < 3; ++i )
@@ -274,8 +279,8 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite
 
     void test_store_and_restore()
     {
-        LUFactorization basis( 3 );
-        LUFactorization otherBasis( 3 );
+        LUFactorization basis( 3, *oracle );
+        LUFactorization otherBasis( 3, *oracle );
 
         double a1[] = { 1, 1, 3 };
         double d1[] = { 0, 0, 0 };
@@ -283,6 +288,14 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite
         TS_ASSERT_THROWS_NOTHING( basis.forwardTransformation( a1, d1 ) );
         basis.pushEtaMatrix( 1, a1 );
 
+        // Save the expected basis after this push
+        double currentBasis[] = {
+            1, 1, 0,
+            0, 1, 0,
+            0, 3, 1
+        };
+        oracle->storeBasis( 3, currentBasis );
+
         // Do a computation using both basis, see that we get the same result.
 
         double a2[] = { 3, 1, 4 };
@@ -327,7 +340,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite
 
     void test_factorization_pivot()
 	{
-		LUFactorization basis( 3 );
+		LUFactorization basis( 3, *oracle );
 		const int nsq = 9;
 
 	    double A[nsq]= {0., 1. , 0.,
@@ -362,7 +375,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite
 	void test_factorization_textbook()
 	{
         // Textbook example
-		LUFactorization basis( 4 );
+		LUFactorization basis( 4, *oracle );
         const int nsq = 16;
         double A[nsq]= {1., 3., -2., 4.,
 						1., 5., -1., 5.,
@@ -438,7 +451,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite
 
 	void test_factorization_as_black_box()
 	{
-		LUFactorization basis( 3 );
+		LUFactorization basis( 3, *oracle );
         const int nsq = 9;
         double A[nsq]= { 1., 2., 4.,
                          4., 5., 7.,
@@ -472,7 +485,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite
 
     void test_factorization_numerical_stability()
     {
-        LUFactorization basis( 3 );
+        LUFactorization basis( 3, *oracle );
 
         double A[] =
             {
@@ -550,13 +563,10 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite
             TS_ASSERT( FloatUtils::areEqual( eta->_column[i], expectedCol3[i] ) );
     }
 
-	void xtest_refactor()
+	void test_refactor()
 	{
-        // TODO: this test fails when the REFACTORIZATION_THRESHOLD is too great (> 10 or so).
-        // Disabling for now.
-
-		LUFactorization basis( 3 );
-		LUFactorization basis2( 3 );
+		LUFactorization basis( 3, *oracle );
+		LUFactorization basis2( 3, *oracle );
 		basis.toggleFactorization( false );
 		int d = 3;
 		unsigned etaCount = GlobalConfiguration::REFACTORIZATION_THRESHOLD + 2;
@@ -574,6 +584,13 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite
                 etaPool[i] = (float)(rand()) / (float)(RAND_MAX);
         }
 
+        double dummyBasis[] = {
+            1, 2, 0,
+            0, 5, 0,
+            0, 8, 9,
+        };
+        oracle->storeBasis( 3, dummyBasis );
+
 		// Generate random etas
 		for ( unsigned i = 0; i < etaCount; ++i )
         {
@@ -586,20 +603,6 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite
 
 		// Check if etas have disappeared
 		TS_ASSERT_EQUALS( basis2.getEtas().size(), etaCount - GlobalConfiguration::REFACTORIZATION_THRESHOLD - 1 );
-		double a[] = {2., -1., 4.};
-		double x1[] = {0., 0., 0.};
-		double y1[] = {0., 0., 0.};
-		double x2[] = {0., 0., 0.};
-		double y2[] = {0., 0., 0.};
-		basis2.forwardTransformation( a, x1 );
-		basis.forwardTransformation( a, y1 );
-		basis2.backwardTransformation( a, x2 );
-		basis.backwardTransformation( a, y2 );
-		for ( unsigned i = 0; i < 3; ++i )
-        {
-			TS_ASSERT( FloatUtils::areEqual( x1[i], y1[i], 0.001 ) );
-			TS_ASSERT( FloatUtils::areEqual( x2[i], y2[i], 0.001 ) );
-		}
 	}
 
     void test_matrix_multiply()
@@ -633,7 +636,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite
 
     void test_invert_B0_fail()
     {
-        LUFactorization basis( 3 );
+        LUFactorization basis( 3, *oracle );
 
         double B0[] = { 1, 0, 0,
                         0, 1, 0,
@@ -654,7 +657,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite
 
     void test_invert_B0()
     {
-        LUFactorization basis( 3 );
+        LUFactorization basis( 3, *oracle );
 
         {
             double B0[] = { 1, 0, 0,
@@ -761,7 +764,7 @@ class LUFactorizationTestSuite : public CxxTest::TestSuite
         }
 
         {
-            LUFactorization basis( 4 );
+            LUFactorization basis( 4, *oracle );
 
             double B0[] = { 1, 1, 1, 0,
                             0, 3, 1, 2,
diff --git a/src/engine/Tableau.cpp b/src/engine/Tableau.cpp
index d34cdd566..91ba6988a 100644
--- a/src/engine/Tableau.cpp
+++ b/src/engine/Tableau.cpp
@@ -218,7 +218,7 @@ void Tableau::setDimensions( unsigned m, unsigned n )
     if ( !_basicStatus )
         throw ReluplexError( ReluplexError::ALLOCATION_FAILED, "Tableau::basicStatus" );
 
-    _basisFactorization = BasisFactorizationFactory::createBasisFactorization( _m );
+    _basisFactorization = BasisFactorizationFactory::createBasisFactorization( _m, *this );
     if ( !_basisFactorization )
         throw ReluplexError( ReluplexError::ALLOCATION_FAILED, "Tableau::basisFactorization" );
 
@@ -1089,7 +1089,7 @@ void Tableau::dumpEquations()
 void Tableau::storeState( TableauState &state ) const
 {
     // Set the dimensions
-    state.setDimensions( _m, _n );
+    state.setDimensions( _m, _n, *this );
 
     // Store matrix A
     memcpy( state._A, _A, sizeof(double) * _n * _m );
@@ -1431,7 +1431,7 @@ void Tableau::addRow()
 
     // Allocate a larger basis factorization
     IBasisFactorization *newBasisFactorization =
-        BasisFactorizationFactory::createBasisFactorization( newM );
+        BasisFactorizationFactory::createBasisFactorization( newM, *this );
     if ( !newBasisFactorization )
         throw ReluplexError( ReluplexError::ALLOCATION_FAILED, "Tableau::newBasisFactorization" );
     delete _basisFactorization;
@@ -1846,6 +1846,13 @@ void Tableau::registerCostFunctionManager( ICostFunctionManager *costFunctionMan
     _costFunctionManager = costFunctionManager;
 }
 
+const double *Tableau::getColumnOfBasis( unsigned column ) const
+{
+    ASSERT( column < _m );
+    unsigned variable = _basicIndexToVariable[column];
+    return _A + ( variable * _m );
+}
+
 //
 // Local Variables:
 // compile-command: "make -C ../.. "
diff --git a/src/engine/Tableau.h b/src/engine/Tableau.h
index 906d1ab86..127651ebf 100644
--- a/src/engine/Tableau.h
+++ b/src/engine/Tableau.h
@@ -13,19 +13,19 @@
 #ifndef __Tableau_h__
 #define __Tableau_h__
 
+#include "IBasisFactorization.h"
 #include "ITableau.h"
 #include "MString.h"
 #include "Map.h"
 #include "Set.h"
 #include "Statistics.h"
 
-class IBasisFactorization;
 class Equation;
 class ICostFunctionManager;
 class PiecewiseLinearCaseSplit;
 class TableauState;
 
-class Tableau : public ITableau
+class Tableau : public ITableau, public IBasisFactorization::BasisColumnOracle
 {
 public:
     Tableau();
@@ -385,6 +385,8 @@ class Tableau : public ITableau
     Equation *getBasisEquation( unsigned row ) const;
     double *getInverseBasisMatrix() const;
 
+    const double *getColumnOfBasis( unsigned column ) const;
+
 private:
     typedef List<VariableWatcher *> VariableWatchers;
     Map<unsigned, VariableWatchers> _variableToWatchers;
diff --git a/src/engine/TableauState.cpp b/src/engine/TableauState.cpp
index fdeb77b14..a4c378939 100644
--- a/src/engine/TableauState.cpp
+++ b/src/engine/TableauState.cpp
@@ -91,7 +91,7 @@ TableauState::~TableauState()
     }
 }
 
-void TableauState::setDimensions( unsigned m, unsigned n )
+void TableauState::setDimensions( unsigned m, unsigned n, const IBasisFactorization::BasisColumnOracle &oracle )
 {
     _m = m;
     _n = n;
@@ -132,7 +132,7 @@ void TableauState::setDimensions( unsigned m, unsigned n )
     if ( !_variableToIndex )
         throw ReluplexError( ReluplexError::ALLOCATION_FAILED, "TableauState::variableToIndex" );
 
-    _basisFactorization = BasisFactorizationFactory::createBasisFactorization( m );
+    _basisFactorization = BasisFactorizationFactory::createBasisFactorization( m, oracle );
     if ( !_basisFactorization )
         throw ReluplexError( ReluplexError::ALLOCATION_FAILED, "TableauState::basisFactorization" );
 }
diff --git a/src/engine/TableauState.h b/src/engine/TableauState.h
index 114cd4af7..b117518e7 100644
--- a/src/engine/TableauState.h
+++ b/src/engine/TableauState.h
@@ -13,11 +13,10 @@
 #ifndef __TableauState_h__
 #define __TableauState_h__
 
+#include "IBasisFactorization.h"
 #include "ITableau.h"
 #include "Set.h"
 
-class IBasisFactorization;
-
 class TableauState
 {
     /*
@@ -37,7 +36,7 @@ class TableauState
     TableauState();
     ~TableauState();
 
-    void setDimensions( unsigned m, unsigned n );
+    void setDimensions( unsigned m, unsigned n, const IBasisFactorization::BasisColumnOracle &oracle );
 
     /*
       The dimensions of matrix A