diff --git a/src/database.rs b/src/database.rs
index dd72b5a..1c1b0d4 100644
--- a/src/database.rs
+++ b/src/database.rs
@@ -104,6 +104,7 @@ use crate::diag::DiagnosticClass;
 use crate::export;
 use crate::formula::Formula;
 use crate::formula::Label;
+use crate::formula::TypeCode;
 use crate::grammar;
 use crate::grammar::Grammar;
 use crate::grammar::StmtParse;
@@ -781,6 +782,15 @@ impl Database {
         )
     }
 
+    /// Returns the typecode for a given label.
+    #[must_use]
+    pub fn label_typecode(&self, label: Label) -> TypeCode {
+        let sref = self
+            .statement_by_label(label)
+            .expect("Invalid label provided to `label_typecode`.");
+        self.name_result().get_atom(sref.math_at(0).slice)
+    }
+
     /// Export an mmp file for a given statement.
     /// Requires: [`Database::name_pass`], [`Database::scope_pass`]
     pub fn export(&self, stmt: &str) {
diff --git a/src/formula.rs b/src/formula.rs
index 92aedff..3bfb5a2 100644
--- a/src/formula.rs
+++ b/src/formula.rs
@@ -49,10 +49,16 @@ pub type Symbol = Atom;
 /// An atom representing a label (nameck suggests `LAtom` for this)
 pub type Label = Atom;
 
-#[derive(Clone, Default)]
+/// An error occurring during unification
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub enum UnificationError {
+    /// Generic unification failure
+    UnificationFailed,
+}
+
 /// A set of substitutions, mapping variables to a formula
 /// We also could have used `dyn Index<&Label, Output=Box<Formula>>`
-#[derive(Debug)]
+#[derive(Clone, Debug, Default)]
 pub struct Substitutions(HashMap<Label, Formula>);
 
 impl Index<Label> for Substitutions {
@@ -65,6 +71,12 @@ impl Index<Label> for Substitutions {
 }
 
 impl Substitutions {
+    /// Creates a new, empty, set of substitutions
+    #[must_use]
+    pub fn new() -> Self {
+        Self::default()
+    }
+
     /// Augment a substitution with a database reference, to produce a [`SubstitutionsRef`].
     /// The resulting object implements [`Debug`].
     #[must_use]
@@ -111,6 +123,13 @@ impl<'a> IntoIterator for &'a Substitutions {
     }
 }
 
+/// A provider for work variables
+/// Work variables are typically used when a new variable appears in an unification, which cannot be immediately assigned.
+pub trait WorkVariableProvider<E> {
+    /// Provide a new work variable for the given typecode
+    fn new_work_variable(&mut self, typecode: TypeCode) -> Result<Label, E>;
+}
+
 /// A [`Substitutions`] reference in the context of a [`Database`].
 /// This allows the values in the [`Substitutions`] to be resolved,
 #[derive(Copy, Clone)]
@@ -150,6 +169,33 @@ pub struct Formula {
 }
 
 impl Formula {
+    /// Creates a formula with just a variable
+    #[must_use]
+    pub fn from_float(label: Label, typecode: TypeCode) -> Self {
+        let mut tree = Tree::default();
+        let root = tree.add_node(label, &[]);
+        let mut variables = Bitset::new();
+        variables.set_bit(root);
+        Formula {
+            typecode,
+            tree: Arc::new(tree),
+            root,
+            variables,
+        }
+    }
+
+    #[inline]
+    /// Iterates through the labels of a formula, depth-first, pre-order.
+    /// Items are the label, and a boolean indicating whether the current label is a variable or not.
+    #[must_use]
+    pub const fn labels_iter(&self) -> LabelIter<'_> {
+        LabelIter {
+            formula: self,
+            stack: vec![],
+            root: Some(self.root),
+        }
+    }
+
     /// Augment a formula with a database reference, to produce a [`FormulaRef`].
     /// The resulting object implements [`Display`], [`Debug`], and [`IntoIterator`].
     #[must_use]
@@ -216,13 +262,15 @@ impl Formula {
     }
 
     /// Unify this formula with the given formula model
-    /// If successful, this returns the substitutions which needs to be made in
+    /// If successful, the provided `substitutions` are completed
+    /// with the substitutions which needs to be made in
     /// `other` in order to match this formula.
-    #[must_use]
-    pub fn unify(&self, other: &Formula) -> Option<Box<Substitutions>> {
-        let mut substitutions = Substitutions(HashMap::default());
-        self.sub_unify(self.root, other, other.root, &mut substitutions)?;
-        Some(Box::new(substitutions))
+    pub fn unify(
+        &self,
+        other: &Formula,
+        substitutions: &mut Substitutions,
+    ) -> Result<(), UnificationError> {
+        self.sub_unify(self.root, other, other.root, substitutions)
     }
 
     /// Unify a sub-formula
@@ -232,18 +280,22 @@ impl Formula {
         other: &Formula,
         other_node_id: NodeId,
         substitutions: &mut Substitutions,
-    ) -> Option<()> {
+    ) -> Result<(), UnificationError> {
         if other.is_variable(other_node_id) {
             // the model formula is a variable, build or match the substitution
             if let Some(formula) = substitutions.0.get(&other.tree[other_node_id]) {
                 // there already is as substitution for that variable, check equality
-                self.sub_eq(node_id, formula, formula.root).then(|| {})
+                if self.sub_eq(node_id, formula, formula.root) {
+                    Ok(())
+                } else {
+                    Err(UnificationError::UnificationFailed)
+                }
             } else {
                 // store the new substitution and succeed
                 substitutions
                     .0
                     .insert(other.tree[other_node_id], self.sub_formula(node_id));
-                Some(())
+                Ok(())
             }
         } else if self.tree[node_id] == other.tree[other_node_id]
             && self.tree.has_children(node_id) == other.tree.has_children(other_node_id)
@@ -256,10 +308,10 @@ impl Formula {
             {
                 self.sub_unify(s_id, other, o_id, substitutions)?;
             }
-            Some(())
+            Ok(())
         } else {
             // formulas differ, we cannot unify.
-            None
+            Err(UnificationError::UnificationFailed)
         }
     }
 
@@ -325,6 +377,46 @@ impl PartialEq for Formula {
     }
 }
 
+/// An iterator through the labels of a formula.
+/// This iterator sequence is depth-first, postfix (post-order).
+/// It provides the label, and a boolean indicating whether the current label is a variable or not.
+#[derive(Debug)]
+pub struct LabelIter<'a> {
+    formula: &'a Formula,
+    stack: Vec<SiblingIter<'a, Label>>,
+    root: Option<NodeId>,
+}
+
+impl<'a> LabelIter<'a> {
+    #[inline]
+    fn visit_children(&mut self, node_id: NodeId) -> (Label, bool) {
+        self.stack.push(self.formula.tree.children_iter(node_id));
+        (
+            self.formula.tree[node_id],
+            self.formula.is_variable(node_id),
+        )
+    }
+}
+
+impl<'a> Iterator for LabelIter<'a> {
+    type Item = (Label, bool);
+
+    fn next(&mut self) -> Option<Self::Item> {
+        if let Some(node_id) = self.root {
+            self.root = None;
+            return Some(self.visit_children(node_id));
+        }
+        loop {
+            let iter = self.stack.last_mut()?;
+            if let Some(node_id) = iter.next() {
+                return Some(self.visit_children(node_id));
+            }
+            // Last sibling reached, pop and iterate
+            self.stack.pop();
+        }
+    }
+}
+
 /// A [`Formula`] reference in the context of a [`Database`].
 /// This allows the values in the [`Formula`] to be resolved,
 #[derive(Copy, Clone)]
@@ -370,12 +462,7 @@ impl<'a> FormulaRef<'a> {
     /// Computes the typecode of the given node
     /// according to the corresponding statement
     fn compute_typecode_at(&self, node_id: NodeId) -> TypeCode {
-        let sref = self
-            .db
-            .statement_by_label(self.formula.tree[node_id])
-            .expect("Formulas shall only contain valid labels. \
-                It shall always be the case when they are obtained from Grammar::parse_formula or StmtParse::get_formula.");
-        self.db.name_result().get_atom(sref.math_at(0).slice)
+        self.db.label_typecode(self.formula.tree[node_id])
     }
 
     /// Convert this formula into an s-expression string.
@@ -406,6 +493,38 @@ impl<'a> FormulaRef<'a> {
         }
     }
 
+    /// Handles the variables present in the formula but not in the substitution list
+    /// The function `f` provided can modify on the fly the substitution list, adding any missing one.
+    pub fn complete_substitutions<E>(
+        &self,
+        substitutions: &mut Substitutions,
+        wvp: &mut impl WorkVariableProvider<E>,
+    ) -> Result<(), E> {
+        self.sub_complete_substitutions(self.formula.root, substitutions, wvp)
+    }
+
+    /// Handles the variables present in the sub-formula but not in the substitution list
+    fn sub_complete_substitutions<E>(
+        &self,
+        node_id: NodeId,
+        substitutions: &mut Substitutions,
+        wvp: &mut impl WorkVariableProvider<E>,
+    ) -> Result<(), E> {
+        if self.is_variable(node_id) {
+            let label = &self.tree[node_id];
+            if substitutions.0.get(label).is_none() {
+                let typecode = self.db.label_typecode(*label);
+                let work_var = wvp.new_work_variable(typecode)?;
+                substitutions.insert(*label, Formula::from_float(work_var, typecode));
+            }
+        } else {
+            for child_node_id in self.tree.children_iter(node_id) {
+                self.sub_complete_substitutions(child_node_id, substitutions, wvp)?;
+            }
+        }
+        Ok(())
+    }
+
     /// Appends this formula to the provided stack buffer.
     ///
     /// The [`ProofBuilder`] structure uses a dense representation of formulas as byte strings,
diff --git a/src/formula_tests.rs b/src/formula_tests.rs
index 0544db3..4f492b5 100644
--- a/src/formula_tests.rs
+++ b/src/formula_tests.rs
@@ -1,4 +1,7 @@
-use crate::grammar_tests::mkdb;
+use crate::{
+    formula::{Substitutions, UnificationError},
+    grammar_tests::mkdb,
+};
 
 const FORMULA_DB: &[u8] = b"
     $c |- wff class ( ) + = 1 2 $.
@@ -32,7 +35,8 @@ fn test_unify() {
     let axiom = stmt_parse
         .get_formula(&db.statement(b"ax-com").unwrap())
         .unwrap();
-    let subst = goal.unify(axiom).unwrap();
+    let mut subst = Substitutions::new();
+    assert!(goal.unify(axiom, &mut subst).is_ok());
     let a = names.lookup_label(b"cA").unwrap().atom;
     let b = names.lookup_label(b"cB").unwrap().atom;
     assert_eq!(subst[a].as_ref(&db).as_sexpr(), "c1");
@@ -50,7 +54,11 @@ fn test_unify_fail() {
     let axiom = stmt_parse
         .get_formula(&db.statement(b"addeq1").unwrap())
         .unwrap();
-    assert!(goal.unify(axiom).is_none());
+    let mut subst = Substitutions::new();
+    assert_eq!(
+        goal.unify(axiom, &mut subst),
+        Err(UnificationError::UnificationFailed)
+    );
 }
 
 #[test]
@@ -66,7 +74,8 @@ fn test_substitute() {
     let axiom = stmt_parse
         .get_formula(&db.statement(b"addeq1.1").unwrap())
         .unwrap();
-    let subst = goal.unify(axiom).unwrap();
+    let mut subst = Substitutions::default();
+    assert!(goal.unify(axiom, &mut subst).is_ok());
     let a = names.lookup_label(b"cA").unwrap().atom;
     let b = names.lookup_label(b"cB").unwrap().atom;
     assert_eq!(subst[a].as_ref(&db).as_sexpr(), "(cadd c1 c2)");
diff --git a/src/grammar.rs b/src/grammar.rs
index 19b0496..22a2ac1 100644
--- a/src/grammar.rs
+++ b/src/grammar.rs
@@ -1191,6 +1191,7 @@ impl Grammar {
         &self,
         symbol_iter: &mut impl Iterator<Item = FormulaToken>,
         expected_typecodes: &[TypeCode],
+        convert_to_provable: bool,
         nset: &Nameset,
     ) -> Result<Formula, StmtParseError> {
         struct StackElement {
@@ -1210,7 +1211,10 @@ impl Grammar {
         let mut stack = vec![];
         loop {
             match *self.nodes.get(e.node_id) {
-                GrammarNode::Leaf { reduce, typecode } => {
+                GrammarNode::Leaf {
+                    reduce,
+                    mut typecode,
+                } => {
                     // We found a leaf: REDUCE
                     Self::do_reduce(&mut formula_builder, reduce, nset);
 
@@ -1243,6 +1247,9 @@ impl Grammar {
                             }
                         } else if symbol_enum.peek().is_none() {
                             // We popped the last element from the stack and we are at the end of the math string, success
+                            if typecode == self.logic_type && convert_to_provable {
+                                typecode = self.provable_type;
+                            }
                             return Ok(formula_builder.build(typecode));
                         } else {
                             // There are still symbols to parse, continue from root
@@ -1263,7 +1270,13 @@ impl Grammar {
                                 {
                                     let reduce = Reduce::new(*label, 1);
                                     Self::do_reduce(&mut formula_builder, reduce, nset);
-                                    return Ok(formula_builder.build(*to_typecode));
+                                    let typecode =
+                                        if *to_typecode == self.logic_type && convert_to_provable {
+                                            self.provable_type
+                                        } else {
+                                            *to_typecode
+                                        };
+                                    return Ok(formula_builder.build(typecode));
                                 }
                             }
                         }
@@ -1408,7 +1421,13 @@ impl Grammar {
         } else {
             typecode.symbol
         };
-        self.parse_formula(&mut symbols, &[expected_typecode], nset)
+        let convert_to_provable = typecode.symbol == self.provable_type;
+        self.parse_formula(
+            &mut symbols,
+            &[expected_typecode],
+            convert_to_provable,
+            nset,
+        )
     }
 
     fn parse_statement(
@@ -1457,8 +1476,13 @@ impl Grammar {
                 }
             })
             .collect();
-        let formula =
-            self.parse_formula(&mut math_string?.into_iter(), &[expected_typecode], nset)?;
+        let convert_to_provable = typecode == self.provable_type;
+        let formula = self.parse_formula(
+            &mut math_string?.into_iter(),
+            &[expected_typecode],
+            convert_to_provable,
+            nset,
+        )?;
         Ok(Some(formula))
     }
 
diff --git a/src/grammar_tests.rs b/src/grammar_tests.rs
index 95c8983..adbbc38 100644
--- a/src/grammar_tests.rs
+++ b/src/grammar_tests.rs
@@ -112,7 +112,12 @@ fn test_parse_formula() {
     let close_parens = make_tok(b")", &names);
     let fmla_vec = vec![a, eq, open_parens, b, plus, a, close_parens];
     let formula = grammar
-        .parse_formula(&mut fmla_vec.clone().into_iter(), &[wff, class], &names)
+        .parse_formula(
+            &mut fmla_vec.clone().into_iter(),
+            &[wff, class],
+            false,
+            &names,
+        )
         .unwrap();
     // Accessing formula using paths to labels
     assert_eq!(
@@ -142,6 +147,14 @@ fn test_parse_formula() {
         .as_ref(&db)
         .iter()
         .eq(fmla_vec.into_iter().map(|t| t.symbol)));
+    // Accessing formula labels
+    let map_label = |(label, is_var)| (as_str(names.atom_name(label)), is_var);
+    let mut labels = formula.labels_iter();
+    assert_eq!(labels.next().map(map_label), Some(("weq", false)));
+    assert_eq!(labels.next().map(map_label), Some(("cA", true)));
+    assert_eq!(labels.next().map(map_label), Some(("cadd", false)));
+    assert_eq!(labels.next().map(map_label), Some(("cB", true)));
+    assert_eq!(labels.next().map(map_label), Some(("cA", true)));
 }
 
 #[test]
@@ -199,7 +212,12 @@ fn test_setvar_as_class() {
     };
     {
         let formula = grammar
-            .parse_formula(&mut vec![x_symbol].into_iter(), &[class_symbol], &names)
+            .parse_formula(
+                &mut vec![x_symbol].into_iter(),
+                &[class_symbol],
+                false,
+                &names,
+            )
             .unwrap();
         assert_eq!(formula.as_ref(&db).as_sexpr(), "(cv vx)");
     }
diff --git a/src/nameck.rs b/src/nameck.rs
index f288b69..836082e 100644
--- a/src/nameck.rs
+++ b/src/nameck.rs
@@ -62,7 +62,7 @@ use std::u32;
 /// a new type like `LAtom` for labels, with `Atom` for symbols only.
 ///
 /// [INC]: /~https://github.com/sorear/smetamath-rs/issues/11
-#[derive(Copy, Clone, Debug, PartialEq, Eq, Default, Hash)]
+#[derive(Copy, Clone, Debug, PartialOrd, Ord, PartialEq, Eq, Default, Hash)]
 pub struct Atom(u32);
 
 // currently we use Vecs for a lot of things in the index.  we might consider
@@ -363,6 +363,47 @@ impl Nameset {
         })
     }
 
+    /// Looks up the float declaration for a math symbol.
+    #[must_use]
+    pub fn lookup_float(&self, symbol: TokenPtr<'_>) -> Option<LookupFloatOwned> {
+        if let Some(syminfo) = self.symbols.get(symbol) {
+            syminfo
+                .float
+                .first()
+                .map(|(addr, (label, _, tcatom))| LookupFloatOwned {
+                    statement_atom: self.get_atom(label),
+                    address: *addr,
+                    typecode_atom: *tcatom,
+                })
+        } else {
+            None
+        }
+    }
+
+    /// Looks up the label declaration for an atom representing a label.
+    /// Atoms are only created and provided for existing labels, so this always returns a value
+    #[must_use]
+    pub fn lookup_label_by_atom(&self, label_atom: Atom) -> LookupLabel {
+        self.lookup_label(self.atom_name(label_atom))
+            .expect("Unknown label atom")
+    }
+
+    /// Looks up up the address and type for an atom representing a math symbol.
+    /// Atoms are only created and provided for existing symbol, so this always returns a value
+    #[must_use]
+    pub fn lookup_symbol_by_atom(&self, symbol_atom: Atom) -> LookupSymbol {
+        self.lookup_symbol(self.atom_name(symbol_atom))
+            .expect("Unknown symbol atom")
+    }
+
+    /// Looks up the float declaration for an atom representing a math symbol.
+    /// Atoms are only created and provided for existing symbol, so this always returns a value
+    #[must_use]
+    pub fn lookup_float_by_atom(&self, symbol_atom: Atom) -> LookupFloatOwned {
+        self.lookup_float(self.atom_name(symbol_atom))
+            .expect("Unknown symbol atom")
+    }
+
     /// Looks up the atom from a $f statement.
     #[must_use]
     pub fn var_atom(&self, sref: StatementRef<'_>) -> Option<Atom> {
@@ -426,6 +467,18 @@ pub struct LookupSymbol {
     pub const_address: Option<TokenAddress>,
 }
 
+/// A representation of that which is stored for each variable in relation to
+/// global $f statements.  All fields apply to the topmost such statement.
+#[derive(Clone, Copy, Debug)]
+pub struct LookupFloatOwned {
+    /// Atom generated for the float statement
+    pub statement_atom: Atom,
+    /// Address of the $f statement.
+    pub address: StatementAddress,
+    /// Atom generated for the typecode.
+    pub typecode_atom: Atom,
+}
+
 /// A representation of that which is stored for each variable in relation to
 /// global $f statements.  All fields apply to the topmost such statement.
 #[derive(Debug)]