16 files changed, 768 insertions, 943 deletions

diff --git a/pyecsca/ec/coordinates.py b/pyecsca/ec/coordinates.py
index 49452c7..8692809 100644
--- a/pyecsca/ec/coordinates.py
+++ b/pyecsca/ec/coordinates.py
@@ -7,8 +7,8 @@ from typing import List, Any, MutableMapping
 
 from public import public
 
-from .formula import (
-    Formula,
+from .formula import Formula
+from .formula.efd import (
     EFDFormula,
     AdditionEFDFormula,
     DoublingEFDFormula,
@@ -55,7 +55,9 @@ class CoordinateModel:
         return f"{self.curve_model.shortname}/{self.name}"
 
     def __repr__(self):
-        return f"{self.__class__.__name__}(\"{self.name}\", curve_model={self.curve_model})"
+        return (
+            f'{self.__class__.__name__}("{self.name}", curve_model={self.curve_model})'
+        )
 
     def __getstate__(self):
         state = self.__dict__.copy()
diff --git a/pyecsca/ec/formula/__init__.py b/pyecsca/ec/formula/__init__.py
new file mode 100644
index 0000000..b6efd8a
--- /dev/null
+++ b/pyecsca/ec/formula/__init__.py
@@ -0,0 +1,4 @@
+""""""
+
+from .base import *
+from .efd import *
diff --git a/pyecsca/ec/formula.py b/pyecsca/ec/formula/base.py
index fff3210..733a82d 100644
--- a/pyecsca/ec/formula.py
+++ b/pyecsca/ec/formula/base.py
@@ -1,24 +1,22 @@
-"""Provides an abstract base class of a formula along with concrete instantiations."""
+"""Provides an abstract base class of a formula."""
 from abc import ABC, abstractmethod
-from ast import parse, Expression
+from ast import Expression
 from functools import cached_property
 
 from astunparse import unparse
-from itertools import product
 from typing import List, Set, Any, ClassVar, MutableMapping, Tuple, Union, Dict
 
-from importlib_resources.abc import Traversable
 from public import public
 from sympy import FF, symbols, Poly, Rational, simplify
-from ..misc.cache import sympify
 
-from .context import ResultAction
-from . import context
-from .error import UnsatisfiedAssumptionError, raise_unsatisified_assumption
-from .mod import Mod, SymbolicMod
-from .op import CodeOp, OpType
-from ..misc.cfg import getconfig
-from ..misc.utils import peval
+from ..context import ResultAction
+from .. import context
+from ..error import UnsatisfiedAssumptionError, raise_unsatisified_assumption
+from ..mod import Mod, SymbolicMod
+from ..op import CodeOp, OpType
+from ...misc.cfg import getconfig
+from ...misc.utils import peval
+from ...misc.cache import sympify
 
 
 @public
@@ -235,7 +233,7 @@ class Formula(ABC):
         :param params: Parameters of the curve.
         :return: The resulting point(s).
         """
-        from .point import Point
+        from ..point import Point
 
         self.__validate_params(field, params)
         self.__validate_points(field, points, params)
@@ -351,93 +349,6 @@ class Formula(ABC):
         )
 
 
-class EFDFormula(Formula):
-    """Formula from the [EFD]_."""
-
-    def __init__(
-        self,
-        meta_path: Traversable,
-        op3_path: Traversable,
-        name: str,
-        coordinate_model: Any,
-    ):
-        self.name = name
-        self.coordinate_model = coordinate_model
-        self.meta = {}
-        self.parameters = []
-        self.assumptions = []
-        self.code = []
-        self.unified = False
-        self.__read_meta_file(meta_path)
-        self.__read_op3_file(op3_path)
-
-    def __read_meta_file(self, path: Traversable):
-        with path.open("rb") as f:
-            line = f.readline().decode("ascii").rstrip()
-            while line:
-                if line.startswith("source"):
-                    self.meta["source"] = line[7:]
-                elif line.startswith("parameter"):
-                    self.parameters.append(line[10:])
-                elif line.startswith("assume"):
-                    self.assumptions.append(
-                        parse(
-                            line[7:].replace("=", "==").replace("^", "**"), mode="eval"
-                        )
-                    )
-                elif line.startswith("unified"):
-                    self.unified = True
-                line = f.readline().decode("ascii").rstrip()
-
-    def __read_op3_file(self, path: Traversable):
-        with path.open("rb") as f:
-            for line in f.readlines():
-                code_module = parse(
-                    line.decode("ascii").replace("^", "**"), str(path), mode="exec"
-                )
-                self.code.append(CodeOp(code_module))
-
-    def __str__(self):
-        return f"{self.coordinate_model!s}/{self.name}"
-
-    @cached_property
-    def input_index(self):
-        return 1
-
-    @cached_property
-    def output_index(self):
-        return max(self.num_inputs + 1, 3)
-
-    @cached_property
-    def inputs(self):
-        return {
-            var + str(i)
-            for var, i in product(
-                self.coordinate_model.variables, range(1, 1 + self.num_inputs)
-            )
-        }
-
-    @cached_property
-    def outputs(self):
-        return {
-            var + str(i)
-            for var, i in product(
-                self.coordinate_model.variables,
-                range(self.output_index, self.output_index + self.num_outputs),
-            )
-        }
-
-    def __eq__(self, other):
-        if not isinstance(other, EFDFormula):
-            return False
-        return (
-            self.name == other.name and self.coordinate_model == other.coordinate_model
-        )
-
-    def __hash__(self):
-        return hash((self.coordinate_model, self.name))
-
-
 @public
 class AdditionFormula(Formula, ABC):
     """Formula that adds two points."""
@@ -448,11 +359,6 @@ class AdditionFormula(Formula, ABC):
 
 
 @public
-class AdditionEFDFormula(AdditionFormula, EFDFormula):
-    pass
-
-
-@public
 class DoublingFormula(Formula, ABC):
     """Formula that doubles a point."""
 
@@ -462,11 +368,6 @@ class DoublingFormula(Formula, ABC):
 
 
 @public
-class DoublingEFDFormula(DoublingFormula, EFDFormula):
-    pass
-
-
-@public
 class TriplingFormula(Formula, ABC):
     """Formula that triples a point."""
 
@@ -476,11 +377,6 @@ class TriplingFormula(Formula, ABC):
 
 
 @public
-class TriplingEFDFormula(TriplingFormula, EFDFormula):
-    pass
-
-
-@public
 class NegationFormula(Formula, ABC):
     """Formula that negates a point."""
 
@@ -490,11 +386,6 @@ class NegationFormula(Formula, ABC):
 
 
 @public
-class NegationEFDFormula(NegationFormula, EFDFormula):
-    pass
-
-
-@public
 class ScalingFormula(Formula, ABC):
     """Formula that somehow scales the point (to a given representative of a projective class)."""
 
@@ -504,11 +395,6 @@ class ScalingFormula(Formula, ABC):
 
 
 @public
-class ScalingEFDFormula(ScalingFormula, EFDFormula):
-    pass
-
-
-@public
 class DifferentialAdditionFormula(Formula, ABC):
     """
     Differential addition formula that adds two points with a known difference.
@@ -522,11 +408,6 @@ class DifferentialAdditionFormula(Formula, ABC):
 
 
 @public
-class DifferentialAdditionEFDFormula(DifferentialAdditionFormula, EFDFormula):
-    pass
-
-
-@public
 class LadderFormula(Formula, ABC):
     """
     Ladder formula for simultaneous addition of two points and doubling of the one of them, with a known difference.
@@ -539,8 +420,3 @@ class LadderFormula(Formula, ABC):
     shortname = "ladd"
     num_inputs = 3
     num_outputs = 2
-
-
-@public
-class LadderEFDFormula(LadderFormula, EFDFormula):
-    pass
diff --git a/pyecsca/ec/formula/efd.py b/pyecsca/ec/formula/efd.py
new file mode 100644
index 0000000..0156fb3
--- /dev/null
+++ b/pyecsca/ec/formula/efd.py
@@ -0,0 +1,142 @@
+""""""
+from functools import cached_property
+from itertools import product
+
+from public import public
+
+from importlib_resources.abc import Traversable
+from typing import Any
+from .base import (
+    Formula,
+    CodeOp,
+    AdditionFormula,
+    DoublingFormula,
+    TriplingFormula,
+    NegationFormula,
+    ScalingFormula,
+    DifferentialAdditionFormula,
+    LadderFormula,
+)
+from ast import parse
+
+
+class EFDFormula(Formula):
+    """Formula from the [EFD]_."""
+
+    def __init__(
+        self,
+        meta_path: Traversable,
+        op3_path: Traversable,
+        name: str,
+        coordinate_model: Any,
+    ):
+        self.name = name
+        self.coordinate_model = coordinate_model
+        self.meta = {}
+        self.parameters = []
+        self.assumptions = []
+        self.code = []
+        self.unified = False
+        self.__read_meta_file(meta_path)
+        self.__read_op3_file(op3_path)
+
+    def __read_meta_file(self, path: Traversable):
+        with path.open("rb") as f:
+            line = f.readline().decode("ascii").rstrip()
+            while line:
+                if line.startswith("source"):
+                    self.meta["source"] = line[7:]
+                elif line.startswith("parameter"):
+                    self.parameters.append(line[10:])
+                elif line.startswith("assume"):
+                    self.assumptions.append(
+                        parse(
+                            line[7:].replace("=", "==").replace("^", "**"), mode="eval"
+                        )
+                    )
+                elif line.startswith("unified"):
+                    self.unified = True
+                line = f.readline().decode("ascii").rstrip()
+
+    def __read_op3_file(self, path: Traversable):
+        with path.open("rb") as f:
+            for line in f.readlines():
+                code_module = parse(
+                    line.decode("ascii").replace("^", "**"), str(path), mode="exec"
+                )
+                self.code.append(CodeOp(code_module))
+
+    def __str__(self):
+        return f"{self.coordinate_model!s}/{self.name}"
+
+    @cached_property
+    def input_index(self):
+        return 1
+
+    @cached_property
+    def output_index(self):
+        return max(self.num_inputs + 1, 3)
+
+    @cached_property
+    def inputs(self):
+        return {
+            var + str(i)
+            for var, i in product(
+                self.coordinate_model.variables, range(1, 1 + self.num_inputs)
+            )
+        }
+
+    @cached_property
+    def outputs(self):
+        return {
+            var + str(i)
+            for var, i in product(
+                self.coordinate_model.variables,
+                range(self.output_index, self.output_index + self.num_outputs),
+            )
+        }
+
+    def __eq__(self, other):
+        if not isinstance(other, EFDFormula):
+            return False
+        return (
+            self.name == other.name and self.coordinate_model == other.coordinate_model
+        )
+
+    def __hash__(self):
+        return hash((self.coordinate_model, self.name))
+
+
+@public
+class AdditionEFDFormula(AdditionFormula, EFDFormula):
+    pass
+
+
+@public
+class DoublingEFDFormula(DoublingFormula, EFDFormula):
+    pass
+
+
+@public
+class TriplingEFDFormula(TriplingFormula, EFDFormula):
+    pass
+
+
+@public
+class NegationEFDFormula(NegationFormula, EFDFormula):
+    pass
+
+
+@public
+class ScalingEFDFormula(ScalingFormula, EFDFormula):
+    pass
+
+
+@public
+class DifferentialAdditionEFDFormula(DifferentialAdditionFormula, EFDFormula):
+    pass
+
+
+@public
+class LadderEFDFormula(LadderFormula, EFDFormula):
+    pass
diff --git a/pyecsca/ec/formula/expand.py b/pyecsca/ec/formula/expand.py
new file mode 100644
index 0000000..d17e80c
--- /dev/null
+++ b/pyecsca/ec/formula/expand.py
@@ -0,0 +1,58 @@
+from typing import List
+
+from .efd import EFDFormula
+from .fliparoo import recursive_fliparoo
+from .graph import ModifiedEFDFormula
+from .metrics import ivs_norm
+from .partitions import reduce_all_adds, expand_all_muls, expand_all_nopower2_muls
+from .switch_sign import generate_switched_formulas
+
+
+def reduce_with_similarity(formulas: List[EFDFormula], norm):
+    efd = list(filter(lambda x: not isinstance(x, ModifiedEFDFormula), formulas))
+    reduced_efd = efd
+    similarities = list(map(norm, efd))
+    for formula in formulas:
+        n = norm(formula)
+        if n in similarities:
+            continue
+        similarities.append(n)
+        reduced_efd.append(formula)
+    return reduced_efd
+
+
+def expand_formula_list(formulas: List[EFDFormula]):
+    extended_efd = reduce_with_similarity(formulas, ivs_norm)
+    print(f"Reduced to {len(extended_efd)} formulas")
+
+    fliparood = sum(list(map(recursive_fliparoo, extended_efd)), [])
+    extended_efd.extend(fliparood)
+    print(f"Fliparoo: {len(extended_efd)} formulas")
+    extended_efd = reduce_with_similarity(extended_efd, ivs_norm)
+    print(f"Reduced to {len(extended_efd)} formulas")
+    # list(map(test_formula, extended_efd))
+
+    switch_signs = sum([list(generate_switched_formulas(f)) for f in extended_efd], [])
+    extended_efd.extend(switch_signs)
+    print(f"Switch signs: {len(extended_efd)} formulas")
+    extended_efd = reduce_with_similarity(extended_efd, ivs_norm)
+    print(f"Reduced to {len(extended_efd)} formulas")
+    # list(map(test_formula, extended_efd))
+
+    extended_efd.extend(list(map(reduce_all_adds, extended_efd)))
+    print(f"Compress adds: {len(extended_efd)} formulas")
+    extended_efd = reduce_with_similarity(extended_efd, ivs_norm)
+    print(f"Reduced to {len(extended_efd)} formulas")
+    # list(map(test_formula, extended_efd))
+
+    extended_efd.extend(list(map(expand_all_muls, extended_efd)))
+    print(f"Expand muls: {len(extended_efd)} formulas")
+    extended_efd = reduce_with_similarity(extended_efd, ivs_norm)
+    print(f"Reduced to {len(extended_efd)} formulas")
+    # list(map(test_formula, extended_efd))
+
+    extended_efd.extend(list(map(expand_all_nopower2_muls, extended_efd)))
+    print(f"Expand muls(!=2^):{len(extended_efd)} formulas")
+    extended_efd = reduce_with_similarity(extended_efd, ivs_norm)
+    print(f"Reduced to {len(extended_efd)} formulas")
+    # list(map(test_formula, extended_efd))
diff --git a/pyecsca/ec/formula_gen/fliparoo.py b/pyecsca/ec/formula/fliparoo.py
index 7e40ca0..5601855 100644
--- a/pyecsca/ec/formula_gen/fliparoo.py
+++ b/pyecsca/ec/formula/fliparoo.py
@@ -1,8 +1,9 @@
-from typing import Dict, Set, Iterator, List, Tuple, Type
-from pyecsca.ec.op import OpType
-from pyecsca.ec.formula_gen.formula_graph import EFDFormulaGraph, Node, CodeOpNode, CodeOp, parse 
-from pyecsca.ec.formula import EFDFormula
-from random import randint 
+from typing import Iterator, List, Tuple, Type, Optional
+from ..op import OpType
+from .graph import EFDFormulaGraph, Node, CodeOpNode, CodeOp, parse
+from .efd import EFDFormula
+from random import randint
+
 
 class Fliparoo:
     """
@@ -12,13 +13,13 @@ class Fliparoo:
         - Neither of N1,...,Nk-1 is an output node
     """
 
-    def __init__(self, chain: List[Node], graph: EFDFormulaGraph):
+    def __init__(self, chain: List[CodeOpNode], graph: EFDFormulaGraph):
         self.verify_chain(chain)
         self.nodes = chain
         self.graph = graph
         self.operator = None
 
-    def verify_chain(self, nodes: List[Node]):
+    def verify_chain(self, nodes: List[CodeOpNode]):
         for i, node in enumerate(nodes[:-1]):
             if node.outgoing_nodes != [nodes[i + 1]]:
                 raise BadFliparoo
@@ -72,7 +73,7 @@ class Fliparoo:
 
 
 class MulFliparoo(Fliparoo):
-    def __init__(self, chain: List[Node], graph: EFDFormulaGraph):
+    def __init__(self, chain: List[CodeOpNode], graph: EFDFormulaGraph):
         super().__init__(chain, graph)
         operations = set(node.op.operator for node in self.nodes)
         if len(operations) != 1 or list(operations)[0] != OpType.Mult:
@@ -81,7 +82,7 @@ class MulFliparoo(Fliparoo):
 
 
 class AddSubFliparoo(Fliparoo):
-    def __init__(self, chain: List[Node], graph: EFDFormulaGraph):
+    def __init__(self, chain: List[CodeOpNode], graph: EFDFormulaGraph):
         super().__init__(chain, graph)
         operations = set(node.op.operator for node in self.nodes)
         if not operations.issubset([OpType.Add, OpType.Sub]):
@@ -89,7 +90,7 @@ class AddSubFliparoo(Fliparoo):
 
 
 class AddFliparoo(Fliparoo):
-    def __init__(self, chain: List[Node], graph: EFDFormulaGraph):
+    def __init__(self, chain: List[CodeOpNode], graph: EFDFormulaGraph):
         super().__init__(chain, graph)
         operations = set(node.op.operator for node in self.nodes)
         if len(operations) != 1 or list(operations)[0] != OpType.Add:
@@ -102,7 +103,7 @@ class BadFliparoo(Exception):
 
 
 def find_fliparoos(
-    graph: EFDFormulaGraph, fliparoo_type: Type[Fliparoo] = None
+    graph: EFDFormulaGraph, fliparoo_type: Optional[Type[Fliparoo]] = None
 ) -> List[Fliparoo]:
     """Finds a list of Fliparoos in a graph"""
     fliparoos = []
@@ -132,8 +133,8 @@ def is_subfliparoo(fliparoo: Tuple[Node], longest_fliparoos: List[Fliparoo]) ->
 
 
 def largest_fliparoo(
-    chain: List[Node], graph: EFDFormula, fliparoo_type: Type[Fliparoo] = None
-) -> Fliparoo:
+    chain: List[Node], graph: EFDFormulaGraph, fliparoo_type: Optional[Type[Fliparoo]] = None
+) -> Optional[Fliparoo]:
     """Finds the largest fliparoo in a list of Nodes"""
     for i in range(len(chain) - 1):
         for j in range(len(chain) - 1, i, -1):
@@ -141,7 +142,7 @@ def largest_fliparoo(
             if fliparoo_type:
                 try:
                     fliparoo_type(subchain, graph)
-                except:
+                except BadFliparoo:
                     continue
             try:
                 return MulFliparoo(subchain, graph)
@@ -223,9 +224,10 @@ class DummyNode(Node):
         pass
 
 
-def generate_fliparood_formulas(formula: EFDFormula, rename: bool = True) -> Iterator[EFDFormula]:
-    graph = EFDFormulaGraph()
-    graph.construct_graph(formula, rename)
+def generate_fliparood_formulas(
+    formula: EFDFormula, rename: bool = True
+) -> Iterator[EFDFormula]:
+    graph = EFDFormulaGraph(formula, rename)
     fliparoos = find_fliparoos(graph)
     for fliparoo in fliparoos:
         for flip_graph in generate_fliparood_graphs(fliparoo):
@@ -236,7 +238,7 @@ def generate_fliparood_graphs(fliparoo: Fliparoo) -> Iterator[EFDFormulaGraph]:
     fliparoo = fliparoo.deepcopy()
     last_str = fliparoo.last.result
     disconnect_fliparoo_outputs(fliparoo)
-    
+
     signed_subgraph = extract_fliparoo_signed_inputs(fliparoo)
 
     # Starting with a single list of unconnected signed nodes
@@ -256,7 +258,7 @@ def generate_fliparood_graphs(fliparoo: Fliparoo) -> Iterator[EFDFormulaGraph]:
         final_signed_node = signed_subgraph.nodes[0]
         if final_signed_node.sign != 1:
             continue
-        final_node = final_signed_node.node
+        final_node: CodeOpNode = final_signed_node.node
 
         opstr = f"{last_str} = {final_node.op.left}{final_node.optype.op_str}{final_node.op.right}"
         final_node.op = CodeOp(parse(opstr))
@@ -266,7 +268,6 @@ def generate_fliparood_graphs(fliparoo: Fliparoo) -> Iterator[EFDFormulaGraph]:
 
 
 def extract_fliparoo_signed_inputs(fliparoo: Fliparoo) -> SignedSubGraph:
-
     graph = fliparoo.graph
     signed_inputs = SignedSubGraph([], graph)
     for node in fliparoo:
@@ -293,19 +294,19 @@ def disconnect_fliparoo_outputs(fliparoo: Fliparoo):
 
 
 def reconnect_fliparoo_outputs(graph: EFDFormulaGraph, last_node: Node):
-    dummy = next(filter(lambda x: isinstance(x,DummyNode),graph.nodes))
+    dummy = next(filter(lambda x: isinstance(x, DummyNode), graph.nodes))
     dummy.reconnect_outgoing_nodes(last_node)
     graph.remove_node(dummy)
-    assert not list(filter(lambda x: isinstance(x,DummyNode),graph.nodes))
+    assert not list(filter(lambda x: isinstance(x, DummyNode), graph.nodes))
 
 
 def combine_all_pairs_signed_nodes(
     signed_subgraph: SignedSubGraph, fliparoo: Fliparoo
 ) -> List[SignedSubGraph]:
     signed_subgraphs = []
-    l = signed_subgraph.components
-    for i in range(l):
-        for j in range(i + 1, l):
+    n_components = signed_subgraph.components
+    for i in range(n_components):
+        for j in range(i + 1, n_components):
             csigned_subgraph = signed_subgraph.deepcopy()
             v, w = csigned_subgraph[i], csigned_subgraph[j]
             combine_signed_nodes(csigned_subgraph, v, w, fliparoo)
@@ -319,8 +320,6 @@ def combine_signed_nodes(
     right_signed_node: SignedNode,
     fliparoo: Fliparoo,
 ):
-
-    result = fliparoo.last.result
     left_node, right_node = left_signed_node.node, right_signed_node.node
     sign = 1
     operator = OpType.Mult
@@ -340,7 +339,7 @@ def combine_signed_nodes(
 
     new_node = CodeOpNode.from_str(
         f"Fliparoo{randint(1,100000)}", left_node.result, operator, right_node.result
-    ) # TODO the randint is not ideal
+    )  # TODO the randint is not ideal
     new_node.incoming_nodes = [left_node, right_node]
     left_node.outgoing_nodes.append(new_node)
     right_node.outgoing_nodes.append(new_node)
@@ -350,3 +349,21 @@ def combine_signed_nodes(
     subgraph.remove_node(left_signed_node)
     subgraph.remove_node(right_signed_node)
     subgraph.add_node(new_node)
+
+
+def recursive_fliparoo(formula, depth=2):
+    all_fliparoos = {0: [formula]}
+    counter = 0
+    while depth > counter:
+        prev_level = all_fliparoos[counter]
+        fliparoo_level = []
+        for flipparood_formula in prev_level:
+            rename = not counter  # rename ivs before the first fliparoo
+            for newly_fliparood in generate_fliparood_formulas(
+                flipparood_formula, rename
+            ):
+                fliparoo_level.append(newly_fliparood)
+        counter += 1
+        all_fliparoos[counter] = fliparoo_level
+
+    return sum(all_fliparoos.values(), [])
diff --git a/pyecsca/ec/formula_gen/formula_graph.py b/pyecsca/ec/formula/graph.py
index 75ebccc..3cb0240 100644
--- a/pyecsca/ec/formula_gen/formula_graph.py
+++ b/pyecsca/ec/formula/graph.py
@@ -1,15 +1,15 @@
-from pyecsca.ec.formula import (
+from .efd import (
     EFDFormula,
     DoublingEFDFormula,
     AdditionEFDFormula,
     LadderEFDFormula,
     DifferentialAdditionEFDFormula,
 )
-from pyecsca.ec.op import CodeOp, OpType
+from ..op import CodeOp, OpType
 import matplotlib.pyplot as plt
 import networkx as nx
 from ast import parse
-from typing import Dict, List, Tuple, Set
+from typing import Dict, List, Tuple, Set, Optional, MutableMapping
 from copy import deepcopy
 from abc import ABC, abstractmethod
 
@@ -21,14 +21,52 @@ class Node(ABC):
         self.output_node = False
         self.input_node = False
 
+    @property
     @abstractmethod
     def label(self) -> str:
         pass
 
+    @property
     @abstractmethod
     def result(self) -> str:
         pass
 
+    @property
+    def is_sub(self) -> bool:
+        return False
+
+    @property
+    def is_mul(self) -> bool:
+        return False
+
+    @property
+    def is_add(self) -> bool:
+        return False
+
+    @property
+    def is_id(self) -> bool:
+        return False
+
+    @property
+    def is_sqr(self) -> bool:
+        return False
+
+    @property
+    def is_pow(self) -> bool:
+        return False
+
+    @property
+    def is_inv(self) -> bool:
+        return False
+
+    @property
+    def is_div(self) -> bool:
+        return False
+
+    @property
+    def is_neg(self) -> bool:
+        return False
+
     @abstractmethod
     def __repr__(self) -> str:
         pass
@@ -43,7 +81,6 @@ class Node(ABC):
 
 
 class ConstantNode(Node):
-
     color = "#b41f44"
 
     def __init__(self, i: int):
@@ -56,14 +93,13 @@ class ConstantNode(Node):
 
     @property
     def result(self) -> str:
-        return self.value
+        return str(self.value)
 
     def __repr__(self) -> str:
         return f"Node({self.value})"
 
 
 class CodeOpNode(Node):
-
     color = "#1f78b4"
 
     def __init__(self, op: CodeOp):
@@ -139,7 +175,6 @@ class CodeOpNode(Node):
 
 
 class InputNode(Node):
-
     color = "#b41f44"
 
     def __init__(self, input: str):
@@ -191,46 +226,45 @@ class ModifiedLadderEFDFormula(LadderEFDFormula, ModifiedEFDFormula):
 
 
 class EFDFormulaGraph:
-    def __init__(self):
-        self.nodes: List = None
-        self.input_nodes: Dict = None
-        self.output_names: Set = None
-        self.roots: List = None
+    nodes: List[Node]
+    input_nodes: MutableMapping[str, InputNode]
+    output_names: Set[str]
+    roots: List[Node]
 
-    def construct_graph(self, formula: EFDFormula, rename=True):
+    def __init__(self, formula: EFDFormula, rename=True):
         self._formula = formula  # TODO remove, its here only for to_EFDFormula
         self.output_names = formula.outputs
         self.input_nodes = {v: InputNode(v) for v in formula_input_variables(formula)}
         self.roots = list(self.input_nodes.values())
         self.nodes = self.roots.copy()
-        discovered_nodes = self.input_nodes.copy()
+        discovered_nodes: Dict[str, Node] = self.input_nodes.copy()
         for op in formula.code:
-            node = CodeOpNode(op)
+            code_node = CodeOpNode(op)
             for side in (op.left, op.right):
                 if side is None:
                     continue
                 if isinstance(side, int):
-                    parent_node = ConstantNode(side)
+                    parent_node: Node = ConstantNode(side)
                     self.nodes.append(parent_node)
                     self.roots.append(parent_node)
                 else:
                     parent_node = discovered_nodes[side]
-                parent_node.outgoing_nodes.append(node)
-                node.incoming_nodes.append(parent_node)
-            self.nodes.append(node)
-            discovered_nodes[op.result] = node
+                parent_node.outgoing_nodes.append(code_node)
+                code_node.incoming_nodes.append(parent_node)
+            self.nodes.append(code_node)
+            discovered_nodes[op.result] = code_node
         # flag output nodes
         for output_name in self.output_names:
             discovered_nodes[output_name].output_node = True
 
         # go through the nodes and make sure that every node is root or has parents
         for node in self.nodes:
-            if not node.incoming_nodes and not node in self.roots:
+            if not node.incoming_nodes and node not in self.roots:
                 self.roots.append(node)
         if rename:
             self.reindex()
 
-    def node_index(self, node: CodeOpNode) -> int:
+    def node_index(self, node: Node) -> int:
         return self.nodes.index(node)
 
     def deepcopy(self):
@@ -250,11 +284,11 @@ class EFDFormulaGraph:
             AdditionEFDFormula: ModifiedAdditionEFDFormula,
             DoublingEFDFormula: ModifiedDoublingEFDFormula,
             DifferentialAdditionEFDFormula: ModifiedDifferentialAdditionEFDFormula,
-            LadderEFDFormula: ModifiedEFDFormula,
+            LadderEFDFormula: ModifiedLadderEFDFormula,
         }
-        if not new_formula.__class__ in set(casting.values()):
+        if new_formula.__class__ not in set(casting.values()):
             new_formula.__class__ = casting[new_formula.__class__]
-        return new_formula
+        return new_formula  # type: ignore
 
     def networkx_graph(self) -> nx.DiGraph:
         graph = nx.DiGraph()
@@ -283,10 +317,10 @@ class EFDFormulaGraph:
             level_counter += 1
         # separate into lists
 
-        level_lists = [[] for _ in range(level_counter)]
+        level_lists: List[List[Node]] = [[] for _ in range(level_counter)]
         discovered = []
         for node, l in reversed(levels):
-            if not node in discovered:
+            if node not in discovered:
                 level_lists[l].append(node)
                 discovered.append(node)
         return level_lists
@@ -304,7 +338,7 @@ class EFDFormulaGraph:
                 )
         return positions
 
-    def draw(self, filename: str = None, figsize: Tuple[int, int] = (12, 12)):
+    def draw(self, filename: Optional[str] = None, figsize: Tuple[int, int] = (12, 12)):
         gnx = self.networkx_graph()
         pos = nx.rescale_layout_dict(self.planar_positions())
         plt.figure(figsize=figsize)
@@ -352,10 +386,10 @@ class EFDFormulaGraph:
         self.nodes.append(node)
 
     def reindex(self):
-        results = {}
+        results: Dict[str, str] = {}
         counter = 0
         for node in self.nodes:
-            if node.input_node or isinstance(node, ConstantNode):
+            if not isinstance(node, CodeOpNode):
                 continue
             op = node.op
             result, left, operator, right = (
@@ -387,6 +421,5 @@ class EFDFormulaGraph:
 
 
 def rename_ivs(formula: EFDFormula):
-    graph = EFDFormulaGraph()
-    graph.construct_graph(formula)
+    graph = EFDFormulaGraph(formula)
     return graph.to_EFDFormula()
diff --git a/pyecsca/ec/formula/metrics.py b/pyecsca/ec/formula/metrics.py
new file mode 100644
index 0000000..a0e42eb
--- /dev/null
+++ b/pyecsca/ec/formula/metrics.py
@@ -0,0 +1,100 @@
+from public import public
+from ...sca.re.zvp import unroll_formula
+from .base import Formula
+import warnings
+from typing import Dict
+from operator import itemgetter, attrgetter
+from ..curve import EllipticCurve
+from ..context import DefaultContext, local
+
+
+@public
+def formula_ivs(formula: Formula):
+    one_unroll = unroll_formula(formula)
+    one_results = {}
+    for name, value in one_unroll:
+        if name in formula.outputs:
+            one_results[name] = value
+    one_polys = set(map(itemgetter(1), one_unroll))
+    return one_polys, set(one_results.values())
+
+
+@public
+def ivs_norm(one: Formula):
+    return formula_ivs(one)[0]
+
+
+@public
+def formula_similarity(one: Formula, other: Formula) -> Dict[str, float]:
+    if one.coordinate_model != other.coordinate_model:
+        warnings.warn("Mismatched coordinate model.")
+
+    one_polys, one_result_polys = formula_ivs(one)
+    other_polys, other_result_polys = formula_ivs(other)
+    return {
+        "output": len(one_result_polys.intersection(other_result_polys))
+        / max(len(one_result_polys), len(other_result_polys)),
+        "ivs": len(one_polys.intersection(other_polys))
+        / max(len(one_polys), len(other_polys)),
+    }
+
+
+@public
+def formula_similarity_abs(one: Formula, other: Formula) -> Dict[str, float]:
+    if one.coordinate_model != other.coordinate_model:
+        warnings.warn("Mismatched coordinate model.")
+
+    one_polys, one_result_polys = formula_ivs(one)
+    other_polys, other_result_polys = formula_ivs(other)
+
+    one_polys = set([f if f.LC() > 0 else -f for f in one_polys])
+    other_polys = set([f if f.LC() > 0 else -f for f in other_polys])
+
+    one_result_polys = set([f if f.LC() > 0 else -f for f in one_result_polys])
+    other_result_polys = set([f if f.LC() > 0 else -f for f in other_result_polys])
+    return {
+        "output": len(one_result_polys.intersection(other_result_polys))
+        / max(len(one_result_polys), len(other_result_polys)),
+        "ivs": len(one_polys.intersection(other_polys))
+        / max(len(one_polys), len(other_polys)),
+    }
+
+
+@public
+def formula_similarity_fuzz(
+    one: Formula, other: Formula, curve: EllipticCurve, samples: int = 1000
+) -> Dict[str, float]:
+    if one.coordinate_model != other.coordinate_model:
+        raise ValueError("Mismatched coordinate model.")
+
+    output_matches = 0.0
+    iv_matches = 0.0
+    for _ in range(samples):
+        Paff = curve.affine_random()
+        Qaff = curve.affine_random()
+        Raff = curve.affine_add(Paff, Qaff)
+        P = Paff.to_model(one.coordinate_model, curve)
+        Q = Qaff.to_model(one.coordinate_model, curve)
+        R = Raff.to_model(one.coordinate_model, curve)
+        inputs = (P, Q, R)[: one.num_inputs]
+        with local(DefaultContext()) as ctx:
+            res_one = one(curve.prime, *inputs, **curve.parameters)
+        action_one = ctx.actions.get_by_index([0])
+        ivs_one = set(
+            map(attrgetter("value"), sum(action_one[0].intermediates.values(), []))
+        )
+        with local(DefaultContext()) as ctx:
+            res_other = other(curve.prime, *inputs, **curve.parameters)
+        action_other = ctx.actions.get_by_index([0])
+        ivs_other = set(
+            map(attrgetter("value"), sum(action_other[0].intermediates.values(), []))
+        )
+        iv_matches += len(ivs_one.intersection(ivs_other)) / max(
+            len(ivs_one), len(ivs_other)
+        )
+        one_coords = set(res_one)
+        other_coords = set(res_other)
+        output_matches += len(one_coords.intersection(other_coords)) / max(
+            len(one_coords), len(other_coords)
+        )
+    return {"output": output_matches / samples, "ivs": iv_matches / samples}
diff --git a/pyecsca/ec/formula_gen/partitions.py b/pyecsca/ec/formula/partitions.py
index ae1413d..28af4b3 100644
--- a/pyecsca/ec/formula_gen/partitions.py
+++ b/pyecsca/ec/formula/partitions.py
@@ -1,34 +1,31 @@
-from pyecsca.ec.op import CodeOp
-from typing import Dict, Set, List
+from typing import List, Any, Generator
 from ast import parse
-from pyecsca.ec.op import OpType
-from pyecsca.ec.formula_gen.formula_graph import (
+from ..op import OpType, CodeOp
+from .graph import (
     EFDFormulaGraph,
     CodeOpNode,
     ConstantNode,
     Node,
 )
-from pyecsca.ec.formula_gen.fliparoo import find_fliparoos, AddFliparoo, MulFliparoo
+from .fliparoo import find_fliparoos, AddFliparoo, MulFliparoo
 from copy import deepcopy
-from pyecsca.ec.formula import EFDFormula
+from .efd import EFDFormula
 
 
 def reduce_all_adds(formula: EFDFormula, rename=True) -> EFDFormula:
-    graph = EFDFormulaGraph()
-    graph.construct_graph(formula, rename=rename)
-    fliparoos = find_single_input_add_fliparoos(graph)
-    for fliparoo in fliparoos:
-        reduce_add_fliparoo(fliparoo, copy=False)
+    graph = EFDFormulaGraph(formula, rename=rename)
+    add_fliparoos = find_single_input_add_fliparoos(graph)
+    for add_fliparoo in add_fliparoos:
+        reduce_add_fliparoo(add_fliparoo, copy=False)
     reduce_all_XplusX(graph)
-    fliparoos = find_constant_mul_fliparoos(graph)
-    for fliparoo in fliparoos:
-        reduce_mul_fliparoo(fliparoo, copy=False)
+    mul_fliparoos = find_constant_mul_fliparoos(graph)
+    for mul_fliparoo in mul_fliparoos:
+        reduce_mul_fliparoo(mul_fliparoo, copy=False)
     return graph.to_EFDFormula()
 
 
 def expand_all_muls(formula: EFDFormula, rename=True) -> EFDFormula:
-    graph = EFDFormulaGraph()
-    graph.construct_graph(formula, rename)
+    graph = EFDFormulaGraph(formula, rename)
     enodes = find_expansion_nodes(graph)
     for enode in enodes:
         expand_mul(graph, enode, copy=False)
@@ -36,8 +33,7 @@ def expand_all_muls(formula: EFDFormula, rename=True) -> EFDFormula:
 
 
 def expand_all_nopower2_muls(formula: EFDFormula, rename=True) -> EFDFormula:
-    graph = EFDFormulaGraph()
-    graph.construct_graph(formula, rename)
+    graph = EFDFormulaGraph(formula, rename)
     enodes = find_expansion_nodes(graph, nopower2=True)
     for enode in enodes:
         expand_mul(graph, enode, copy=False)
@@ -75,7 +71,7 @@ def find_constant_mul_fliparoos(graph: EFDFormulaGraph) -> List[MulFliparoo]:
             if len(nonconstant_inputs) != 1:
                 continue
             inode = nonconstant_inputs[0]
-            if not inode in fliparoo.first.incoming_nodes:
+            if inode not in fliparoo.first.incoming_nodes:
                 continue
             if not sum(
                 1
@@ -92,7 +88,7 @@ def find_constant_mul_fliparoos(graph: EFDFormulaGraph) -> List[MulFliparoo]:
 
 
 def find_expansion_nodes(graph: EFDFormulaGraph, nopower2=False) -> List[Node]:
-    expansion_nodes = []
+    expansion_nodes: List[Node] = []
     for node in graph.nodes:
         if not isinstance(node, CodeOpNode) or not node.is_mul:
             continue
@@ -120,7 +116,7 @@ def reduce_all_XplusX(graph: EFDFormulaGraph):
     graph.update()
 
 
-def find_all_XplusX(graph) -> List[Node]:
+def find_all_XplusX(graph) -> List[CodeOpNode]:
     adds = []
     for node in graph.nodes:
         if not isinstance(node, CodeOpNode) or not node.is_add:
@@ -256,8 +252,8 @@ class Partition:
     def __eq__(self, other):
         if self.value != other.value:
             return False
-        if self.final or other.final:
-            return self.final == other.final
+        if self.is_final or other.is_final:
+            return self.is_final == other.is_final
         l, r = self.parts
         lo, ro = other.parts
         return (l == lo and r == ro) or (l == ro and r == lo)
@@ -272,13 +268,14 @@ def compute_partitions(n: int) -> List[Partition]:
                 partitions.append(partition_dp + partition_n_dp)
     # remove duplicates
     result = []
-    [result.append(p) for p in partitions if p not in result]
+    for p in partitions:
+        if p not in result:
+            result.append(p)
     return result
 
 
 def generate_partitioned_formulas(formula: EFDFormula, rename=True):
-    graph = EFDFormulaGraph()
-    graph.construct_graph(formula, rename)
+    graph = EFDFormulaGraph(formula, rename)
     enodes = find_expansion_nodes(graph)
     for enode in enodes:
         for part_graph in generate_all_node_partitions(graph, enode):
@@ -287,8 +284,7 @@ def generate_partitioned_formulas(formula: EFDFormula, rename=True):
 
 def generate_all_node_partitions(
     original_graph: EFDFormulaGraph, node: Node
-) -> EFDFormulaGraph:
-
+) -> Generator[EFDFormulaGraph, Any, None]:
     const_par = next(filter(lambda x: isinstance(x, ConstantNode), node.incoming_nodes))
     const_par_value = const_par.value
 
diff --git a/pyecsca/ec/formula_gen/switch_sign.py b/pyecsca/ec/formula/switch_sign.py
index 4ecf3a1..b8a0121 100644
--- a/pyecsca/ec/formula_gen/switch_sign.py
+++ b/pyecsca/ec/formula/switch_sign.py
@@ -1,17 +1,15 @@
-from pyecsca.ec.op import CodeOp
-from typing import Dict, Set, Iterator, List, Tuple, Type
+from typing import Dict, Iterator, List
 from ast import parse
-from pyecsca.ec.op import OpType
-from pyecsca.ec.formula_gen.formula_graph import EFDFormulaGraph, ConstantNode, Node
+from ..op import OpType, CodeOp
+from .graph import EFDFormulaGraph, ConstantNode, Node, CodeOpNode
 from itertools import chain, combinations
-from pyecsca.ec.formula import EFDFormula
+from .efd import EFDFormula
 
 
 def generate_switched_formulas(
     formula: EFDFormula, rename=True
 ) -> Iterator[EFDFormula]:
-    graph = EFDFormulaGraph()
-    graph.construct_graph(formula, rename)
+    graph = EFDFormulaGraph(formula, rename)
     for node_combination in subnode_lists(graph):
         try:
             yield switch_sign(graph, node_combination).to_EFDFormula()
@@ -61,7 +59,9 @@ class BadSignSwitch(Exception):
     pass
 
 
-def switch_sign_propagate(node: Node, variable: str, output_signs: Dict[Node, str]):
+def switch_sign_propagate(
+    node: CodeOpNode, variable: str, output_signs: Dict[Node, str]
+):
     if node.is_add:
         if variable == node.incoming_nodes[1].result:
             node.op = change_operator(node.op, OpType.Sub)
diff --git a/pyecsca/ec/formula_gen/formula_gen.py b/pyecsca/ec/formula_gen/formula_gen.py
deleted file mode 100644
index 35f452f..0000000
--- a/pyecsca/ec/formula_gen/formula_gen.py
+++ /dev/null
@@ -1,131 +0,0 @@
-from pyecsca.ec.model import ShortWeierstrassModel, MontgomeryModel, TwistedEdwardsModel
-from pyecsca.ec.formula_gen.test import load_efd_formulas, load_library_formulas
-from pyecsca.ec.formula_gen.formula_graph import ModifiedEFDFormula
-from pyecsca.ec.formula_gen.fliparoo import generate_fliparood_formulas
-from pyecsca.ec.formula_gen.switch_sign import generate_switched_formulas
-from pyecsca.ec.formula_gen.partitions import (
-    reduce_all_adds,
-    expand_all_muls,
-    expand_all_nopower2_muls,
-)
-from pyecsca.ec.formula import EFDFormula, Formula
-from typing import List, Dict
-from operator import itemgetter
-from tqdm.notebook import tqdm
-from pyecsca.sca.re.zvp import unroll_formula
-import warnings
-from pyecsca.ec.formula_gen.test import test_formula
-
-
-def main():
-    efd = load_efd_formulas("projective", ShortWeierstrassModel())
-
-    extended_efd = list(efd.values())
-    print(f"{len(extended_efd)} formulas")
-    extended_efd = reduce_with_similarity(extended_efd, iv_set)
-    print(f"Reduced to {len(extended_efd)} formulas")
-
-    fliparood = sum(list(map(recursive_fliparoo, extended_efd)), [])
-    extended_efd.extend(fliparood)
-    print(f"Fliparoo: {len(extended_efd)} formulas")
-    extended_efd = reduce_with_similarity(extended_efd, iv_set)
-    print(f"Reduced to {len(extended_efd)} formulas")
-    list(map(test_formula, extended_efd))
-
-    switch_signs = sum([list(generate_switched_formulas(f)) for f in extended_efd], [])
-    extended_efd.extend(switch_signs)
-    print(f"Switch signs: {len(extended_efd)} formulas")
-    extended_efd = reduce_with_similarity(extended_efd, iv_set)
-    print(f"Reduced to {len(extended_efd)} formulas")
-    list(map(test_formula, extended_efd))
-
-    extended_efd.extend(list(map(reduce_all_adds, extended_efd)))
-    print(f"Compress adds: {len(extended_efd)} formulas")
-    extended_efd = reduce_with_similarity(extended_efd, iv_set)
-    print(f"Reduced to {len(extended_efd)} formulas")
-    list(map(test_formula, extended_efd))
-
-    extended_efd.extend(list(map(expand_all_muls, extended_efd)))
-    print(f"Expand muls: {len(extended_efd)} formulas")
-    extended_efd = reduce_with_similarity(extended_efd, iv_set)
-    print(f"Reduced to {len(extended_efd)} formulas")
-    list(map(test_formula, extended_efd))
-
-    extended_efd.extend(list(map(expand_all_nopower2_muls, extended_efd)))
-    print(f"Expand muls(!=2^):{len(extended_efd)} formulas")
-    extended_efd = reduce_with_similarity(extended_efd, iv_set)
-    print(f"Reduced to {len(extended_efd)} formulas")
-    list(map(test_formula, extended_efd))
-
-    return extended_efd
-
-
-def reduce_with_similarity(formulas: List[EFDFormula], norm):
-    efd = list(filter(lambda x: not isinstance(x, ModifiedEFDFormula), formulas))
-    reduced_efd = efd
-    similarities = list(map(norm, efd))
-    for formula in formulas:
-        n = norm(formula)
-        if n in similarities:
-            continue
-        similarities.append(n)
-        reduced_efd.append(formula)
-    return reduced_efd
-
-
-def formula_similarity(one: Formula, other: Formula) -> Dict[str, float]:
-    if one.coordinate_model != other.coordinate_model:
-        warnings.warn("Mismatched coordinate model.")
-
-    one_unroll = unroll_formula(one)
-    other_unroll = unroll_formula(other)
-    one_results = {}
-    for name, value in one_unroll:
-        if name in one.outputs:
-            one_results[name] = value
-    other_results = {}
-    for name, value in other_unroll:
-        if name in other.outputs:
-            other_results[name] = value
-    one_result_polys = set(one_results.values())
-    other_result_polys = set(other_results.values())
-    one_polys = set(map(itemgetter(1), one_unroll))
-    other_polys = set(map(itemgetter(1), other_unroll))
-    return {
-        "output": len(one_result_polys.intersection(other_result_polys))
-        / max(len(one_result_polys), len(other_result_polys)),
-        "ivs": len(one_polys.intersection(other_polys))
-        / max(len(one_polys), len(other_polys)),
-    }
-
-
-def iv_set(one: Formula):
-    one_unroll = unroll_formula(one)
-    one_results = {}
-    for name, value in one_unroll:
-        if name in one.outputs:
-            one_results[name] = value
-    one_polys = set(map(itemgetter(1), one_unroll))
-    return one_polys
-
-
-def recursive_fliparoo(formula, depth=2):
-    all_fliparoos = {0: [formula]}
-    counter = 0
-    while depth > counter:
-        prev_level = all_fliparoos[counter]
-        fliparoo_level = []
-        for flipparood_formula in prev_level:
-            rename = not counter  # rename ivs before the first fliparoo
-            for newly_fliparood in generate_fliparood_formulas(
-                flipparood_formula, rename
-            ):
-                fliparoo_level.append(newly_fliparood)
-        counter += 1
-        all_fliparoos[counter] = fliparoo_level
-
-    return sum(all_fliparoos.values(), [])
-
-
-if __name__ == "__main__":
-    main()
diff --git a/pyecsca/ec/formula_gen/metrics.py b/pyecsca/ec/formula_gen/metrics.py
deleted file mode 100644
index 45f51c2..0000000
--- a/pyecsca/ec/formula_gen/metrics.py
+++ /dev/null
@@ -1,154 +0,0 @@
-from pyecsca.sca.re.zvp import unroll_formula
-from pyecsca.ec.formula import (
-    EFDFormula,
-    AdditionEFDFormula,
-    Formula,
-    LadderEFDFormula,
-    DoublingEFDFormula,
-)
-import warnings
-from typing import Dict
-from operator import itemgetter, attrgetter
-from pyecsca.ec.curve import EllipticCurve
-from pyecsca.ec.context import DefaultContext, local
-
-
-def formula_similarity(one: Formula, other: Formula) -> Dict[str, float]:
-    if one.coordinate_model != other.coordinate_model:
-        warnings.warn("Mismatched coordinate model.")
-
-    one_unroll = unroll_formula(one)
-    other_unroll = unroll_formula(other)
-    one_results = {}
-    for name, value in one_unroll:
-        if name in one.outputs:
-            one_results[name] = value
-    other_results = {}
-    for name, value in other_unroll:
-        if name in other.outputs:
-            other_results[name] = value
-    one_result_polys = set(one_results.values())
-    other_result_polys = set(other_results.values())
-    one_polys = set(map(itemgetter(1), one_unroll))
-    other_polys = set(map(itemgetter(1), other_unroll))
-    return {
-        "output": len(one_result_polys.intersection(other_result_polys))
-        / max(len(one_result_polys), len(other_result_polys)),
-        "ivs": len(one_polys.intersection(other_polys))
-        / max(len(one_polys), len(other_polys)),
-    }
-
-
-def ivs_norm(one: Formula):
-    one_unroll = unroll_formula(one)
-    one_results = {}
-    for name, value in one_unroll:
-        if name in one.outputs:
-            one_results[name] = value
-    one_polys = set(map(itemgetter(1), one_unroll))
-    return one_polys
-
-
-def formula_similarity_abs(one: Formula, other: Formula) -> Dict[str, float]:
-    if one.coordinate_model != other.coordinate_model:
-        warnings.warn("Mismatched coordinate model.")
-
-    one_unroll = unroll_formula(one)
-    other_unroll = unroll_formula(other)
-    one_results = {}
-    for name, value in one_unroll:
-        if name in one.outputs:
-            one_results[name] = value
-    other_results = {}
-    for name, value in other_unroll:
-        if name in other.outputs:
-            other_results[name] = value
-    one_result_polys = set(one_results.values())
-    other_result_polys = set(other_results.values())
-    one_polys = set(map(itemgetter(1), one_unroll))
-    other_polys = set(map(itemgetter(1), other_unroll))
-
-    one_polys = set([f if f.LC() > 0 else -f for f in one_polys])
-    # one_polys.update(set(-f for f in one_polys))
-    other_polys = set([f if f.LC() > 0 else -f for f in other_polys])
-    # other_polys.update(set(-f for f in other_polys))
-
-    return {
-        "output": len(one_result_polys.intersection(other_result_polys))
-        / max(len(one_result_polys), len(other_result_polys)),
-        "ivs": len(one_polys.intersection(other_polys))
-        / max(len(one_polys), len(other_polys)),
-    }
-
-
-def formula_similarity_fuzz(
-    one: Formula, other: Formula, curve: EllipticCurve, samples: int = 1000
-) -> Dict[str, float]:
-    if one.coordinate_model != other.coordinate_model:
-        raise ValueError("Mismatched coordinate model.")
-
-    output_matches = 0.0
-    iv_matches = 0.0
-    for _ in range(samples):
-        P = curve.affine_random().to_model(one.coordinate_model, curve)
-        Q = curve.affine_random().to_model(other.coordinate_model, curve)
-        with local(DefaultContext()) as ctx:
-            res_one = one(curve.prime, P, Q, **curve.parameters)
-        action_one = ctx.actions.get_by_index([0])
-        ivs_one = set(
-            map(attrgetter("value"), sum(action_one[0].intermediates.values(), []))
-        )
-        with local(DefaultContext()) as ctx:
-            res_other = other(curve.prime, P, Q, **curve.parameters)
-        action_other = ctx.actions.get_by_index([0])
-        ivs_other = set(
-            map(attrgetter("value"), sum(action_other[0].intermediates.values(), []))
-        )
-        iv_matches += len(ivs_one.intersection(ivs_other)) / max(
-            len(ivs_one), len(ivs_other)
-        )
-        one_coords = set(res_one)
-        other_coords = set(res_other)
-        output_matches += len(one_coords.intersection(other_coords)) / max(
-            len(one_coords), len(other_coords)
-        )
-    return {"output": output_matches / samples, "ivs": iv_matches / samples}
-
-
-def formula_similarity_fuzz2(
-    one: Formula, other: Formula, curve: EllipticCurve, samples: int = 1000
-) -> Dict[str, float]:
-    if one.coordinate_model != other.coordinate_model:
-        raise ValueError("Mismatched coordinate model.")
-
-    output_matches = 0.0
-    iv_matches = 0.0
-    for _ in range(samples):
-        Paff = curve.affine_random()
-        Qaff = curve.affine_random()
-        Raff = curve.affine_add(Paff, Qaff)
-        P = Paff.to_model(one.coordinate_model, curve)
-        Q = Qaff.to_model(one.coordinate_model, curve)
-        R = Raff.to_model(one.coordinate_model, curve)
-        inputs = (P, Q, R)[: one.num_inputs]
-        with local(DefaultContext()) as ctx:
-            res_one = one(curve.prime, *inputs, **curve.parameters)
-        action_one = ctx.actions.get_by_index([0])
-        ivs_one = set(
-            map(attrgetter("value"), sum(action_one[0].intermediates.values(), []))
-        )
-        with local(DefaultContext()) as ctx:
-            res_other = other(curve.prime, *inputs, **curve.parameters)
-        action_other = ctx.actions.get_by_index([0])
-        ivs_other = set(
-            map(attrgetter("value"), sum(action_other[0].intermediates.values(), []))
-        )
-        iv_matches += len(ivs_one.intersection(ivs_other)) / max(
-            len(ivs_one), len(ivs_other)
-        )
-        one_coords = set(res_one)
-        other_coords = set(res_other)
-        output_matches += len(one_coords.intersection(other_coords)) / max(
-            len(one_coords), len(other_coords)
-        )
-    return {"output": output_matches / samples, "ivs": iv_matches / samples}
diff --git a/pyecsca/ec/formula_gen/test.py b/pyecsca/ec/formula_gen/test.py
deleted file mode 100644
index 8517909..0000000
--- a/pyecsca/ec/formula_gen/test.py
+++ /dev/null
@@ -1,324 +0,0 @@
-from pyecsca.ec.params import get_params
-from importlib_resources import as_file
-from pyecsca.ec.formula import (
-    AdditionFormula,
-    LadderFormula,
-    DoublingFormula,
-    AdditionEFDFormula,
-    DoublingEFDFormula,
-    LadderEFDFormula,
-)
-
-from pathlib import Path
-from pyecsca.ec.model import ShortWeierstrassModel, MontgomeryModel, TwistedEdwardsModel
-from pyecsca.ec.formula_gen.fliparoo import generate_fliparood_formulas
-from pyecsca.ec.formula_gen.switch_sign import generate_switched_formulas
-from pyecsca.ec.formula_gen.partitions import (
-    reduce_all_adds,
-    expand_all_muls,
-    expand_all_nopower2_muls,
-)
-from pyecsca.ec.formula_gen.formula_graph import rename_ivs
-
-
-def main():
-    for name, lib_formula in load_library_formulas().items():
-        print(name)
-        test_formula_graph(lib_formula, library=True)
-        test_formula(lib_formula, library=True)
-        test_fliparood_formula(lib_formula, library=True)
-        test_switch_sign(lib_formula, library=True)
-        test_reductions(lib_formula, library=True)
-        test_expansions(lib_formula, library=True)
-    for name, formula in load_efd_formulas(
-        "projective", ShortWeierstrassModel()
-    ).items():
-        print(name)
-        test_fliparood_formula(formula)
-        test_switch_sign(formula)
-        test_reductions(formula)
-        test_expansions(formula)
-    print("All good.")
-
-
-def test_formula_graph(formula, library=False):
-    test_formula(rename_ivs(formula), library)
-
-
-def test_switch_sign(formula, library=False):
-    for switch_formula in generate_switched_formulas(formula):
-        test_formula(switch_formula, library)
-
-
-def test_fliparood_formula(formula, library=False):
-    for fliparood in generate_fliparood_formulas(formula):
-        test_formula(fliparood, library)
-
-
-def test_reductions(formula, library=False):
-    test_formula(reduce_all_adds(formula), library)
-
-
-def test_expansions(formula, library=False):
-    test_formula(expand_all_muls(formula), library)
-    test_formula(expand_all_nopower2_muls(formula), library)
-
-
-def test_formula(formula, library=False):
-    try:
-        test_formula0(formula, library)
-    except AssertionError:
-        print(formula.name)
-
-
-def test_formula0(formula, library=False):
-    coordinate_model = formula.coordinate_model
-    param_spec = choose_curve(coordinate_model, formula.name, library)
-    params = get_params(*param_spec, coordinate_model.name)
-    scale = coordinate_model.formulas.get("z", None)
-    if scale is None:
-        scale = coordinate_model.formulas.get("scale", None)
-
-    formula_type = formula.__class__
-    for _ in range(10):
-        Paff = params.curve.affine_random()
-        P2aff = params.curve.affine_double(Paff)
-        Qaff = params.curve.affine_random()
-        Q2aff = params.curve.affine_double(Qaff)
-        Raff = params.curve.affine_add(Paff, Qaff)
-        R2aff = params.curve.affine_double(Raff)
-        QRaff = params.curve.affine_add(Qaff, Raff)
-        P = Paff.to_model(coordinate_model, params.curve)
-        P2 = P2aff.to_model(coordinate_model, params.curve)
-        Q = Qaff.to_model(coordinate_model, params.curve)
-        Q2 = Q2aff.to_model(coordinate_model, params.curve)
-        R = Raff.to_model(coordinate_model, params.curve)
-        R2 = R2aff.to_model(coordinate_model, params.curve)  # noqa
-        QR = QRaff.to_model(coordinate_model, params.curve)
-        inputs = (P, Q, R)[: formula.num_inputs]
-        res = formula(params.curve.prime, *inputs, **params.curve.parameters)
-        if issubclass(formula_type, AdditionFormula):
-            try:
-                assert res[0].to_affine() == Raff
-            except NotImplementedError:
-                assert (
-                    scale(params.curve.prime, res[0], **params.curve.parameters)[0] == R
-                )
-        elif issubclass(formula_type, DoublingFormula):
-            try:
-                assert res[0].to_affine() == P2aff
-            except NotImplementedError:
-                assert (
-                    scale(params.curve.prime, res[0], **params.curve.parameters)[0]
-                    == P2
-                )
-        elif issubclass(formula_type, LadderFormula):
-            try:
-                assert res[0].to_affine() == Q2aff
-                assert res[1].to_affine() == QRaff
-            except NotImplementedError:
-                # print(scale(params.curve.prime, res[0], **params.curve.parameters)[0])
-                # print(scale(params.curve.prime, res[1], **params.curve.parameters)[0])
-                # print(P)
-                # print(Q)
-                # print(R)
-                # print(P2)
-                # print(Q2)
-                # print(R2)
-                # print(QR)
-                # print("------------------------------------")
-                assert (
-                    scale(params.curve.prime, res[1], **params.curve.parameters)[0]
-                    == QR
-                )
-                assert (
-                    scale(params.curve.prime, res[0], **params.curve.parameters)[0]
-                    == Q2
-                )
-
-
-def load_efd_formulas(coordinate_name, model):
-    formulas = model.coordinates[coordinate_name].formulas
-    return {name: f for name, f in formulas.items() if "add" in name or "dbl" in name}
-
-
-def load_library_formulas(coordinates=None):
-    libs_dict = {}
-    for name, model, coords, _, formula_type in LIBRARY_FORMULAS:
-        if coordinates is not None and coordinates != coords:
-            continue
-        coordinate_model = model().coordinates[coords]
-        lib_path = Path("test/data/formulas")  # Path("../../../test/data/formulas")
-        with as_file(lib_path.joinpath(name)) as meta_path, as_file(
-            lib_path.joinpath(name + ".op3")
-        ) as op3_path:
-            libs_dict[name] = formula_type(meta_path, op3_path, name, coordinate_model)
-    return libs_dict
-
-
-LIBRARY_FORMULAS = [
-    [
-        "add-bc-r1rv76-jac",
-        ShortWeierstrassModel,
-        "jacobian",
-        ("secg", "secp128r1"),
-        AdditionEFDFormula,
-    ],
-    [
-        "add-bc-r1rv76-mod",
-        ShortWeierstrassModel,
-        "modified",
-        ("secg", "secp128r1"),
-        AdditionEFDFormula,
-    ],
-    [
-        "dbl-bc-r1rv76-jac",
-        ShortWeierstrassModel,
-        "jacobian",
-        ("secg", "secp128r1"),
-        DoublingEFDFormula,
-    ],
-    [
-        "dbl-bc-r1rv76-mod",
-        ShortWeierstrassModel,
-        "modified",
-        ("secg", "secp128r1"),
-        DoublingEFDFormula,
-    ],
-    [
-        "dbl-bc-r1rv76-x25519",
-        MontgomeryModel,
-        "xz",
-        ("other", "Curve25519"),
-        DoublingEFDFormula,
-    ],
-    [
-        "ladd-bc-r1rv76-x25519",
-        MontgomeryModel,
-        "xz",
-        ("other", "Curve25519"),
-        LadderEFDFormula,
-    ],
-    [
-        "dbl-boringssl-p224",
-        ShortWeierstrassModel,
-        "jacobian-3",
-        ("secg", "secp224r1"),
-        DoublingEFDFormula,
-    ],
-    [
-        "add-boringssl-p224",
-        ShortWeierstrassModel,
-        "jacobian-3",
-        ("secg", "secp224r1"),
-        AdditionEFDFormula,
-    ],
-    [
-        "add-libressl-v382",
-        ShortWeierstrassModel,
-        "jacobian",
-        ("secg", "secp128r1"),
-        AdditionEFDFormula,
-    ],
-    [
-        "dbl-libressl-v382",
-        ShortWeierstrassModel,
-        "jacobian",
-        ("secg", "secp128r1"),
-        DoublingEFDFormula,
-    ],
-    [
-        "dbl-secp256k1-v040",
-        ShortWeierstrassModel,
-        "jacobian",
-        ("secg", "secp256k1"),
-        DoublingEFDFormula,
-    ],
-    [
-        "add-openssl-z256",
-        ShortWeierstrassModel,
-        "jacobian-3",
-        ("secg", "secp256r1"),
-        AdditionEFDFormula,
-    ],
-    [
-        "add-openssl-z256a",
-        ShortWeierstrassModel,
-        "jacobian-3",
-        ("secg", "secp256r1"),
-        AdditionEFDFormula,
-    ],
-    [
-        "ladd-openssl-x25519",
-        MontgomeryModel,
-        "xz",
-        ("other", "Curve25519"),
-        LadderEFDFormula,
-    ],
-    [
-        "ladd-hacl-x25519",
-        MontgomeryModel,
-        "xz",
-        ("other", "Curve25519"),
-        LadderEFDFormula,
-    ],
-    [
-        "dbl-hacl-x25519",
-        MontgomeryModel,
-        "xz",
-        ("other", "Curve25519"),
-        DoublingEFDFormula,
-    ],
-    [
-        "dbl-sunec-v21",
-        ShortWeierstrassModel,
-        "projective-3",
-        ("secg", "secp256r1"),
-        DoublingEFDFormula,
-    ],
-    [
-        "add-sunec-v21",
-        ShortWeierstrassModel,
-        "projective-3",
-        ("secg", "secp256r1"),
-        AdditionEFDFormula,
-    ],
-    [
-        "add-sunec-v21-ed25519",
-        TwistedEdwardsModel,
-        "extended",
-        ("other", "Ed25519"),
-        AdditionEFDFormula,
-    ],
-    [
-        "dbl-sunec-v21-ed25519",
-        TwistedEdwardsModel,
-        "extended",
-        ("other", "Ed25519"),
-        DoublingEFDFormula,
-    ],
-    [
-        "ladd-rfc7748",
-        MontgomeryModel,
-        "xz",
-        ("other", "Curve25519"),
-        LadderEFDFormula,
-    ],
-]
-
-
-def choose_curve(coordinate_model, name, library):
-    if library:
-        return next(filter(lambda x: x[0] == name, LIBRARY_FORMULAS))[3]
-    model = coordinate_model.curve_model
-    if model.__class__ == ShortWeierstrassModel:
-        return ("secg", "secp128r1")
-    if model.__class__ == MontgomeryModel:
-        return ("other", "Curve25519")
-    if model.__class__ == TwistedEdwardsModel:
-        return ("other", "Ed25519")
-    raise NotImplementedError(model)
-
-
-if __name__ == "__main__":
-    main()
diff --git a/pyecsca/sca/re/structural.py b/pyecsca/sca/re/structural.py
index c79e604..f3b0d32 100644
--- a/pyecsca/sca/re/structural.py
+++ b/pyecsca/sca/re/structural.py
@@ -1,77 +1 @@
 """"""
-import warnings
-from typing import Dict
-from public import public
-
-from ...ec.curve import EllipticCurve
-from ...ec.formula import Formula
-from ...ec.context import DefaultContext, local
-from .zvp import unroll_formula
-from operator import itemgetter, attrgetter
-
-
-@public
-def formula_similarity(one: Formula, other: Formula) -> Dict[str, float]:
-    if one.coordinate_model != other.coordinate_model:
-        warnings.warn("Mismatched coordinate model.")
-
-    one_unroll = unroll_formula(one)
-    other_unroll = unroll_formula(other)
-    one_results = {}
-    for name, value in one_unroll:
-        if name in one.outputs:
-            one_results[name] = value
-    other_results = {}
-    for name, value in other_unroll:
-        if name in other.outputs:
-            other_results[name] = value
-    one_result_polys = set(one_results.values())
-    other_result_polys = set(other_results.values())
-    one_polys = set(map(itemgetter(1), one_unroll))
-    other_polys = set(map(itemgetter(1), other_unroll))
-    return {
-        "output": len(one_result_polys.intersection(other_result_polys))
-        / max(len(one_result_polys), len(other_result_polys)),
-        "ivs": len(one_polys.intersection(other_polys))
-        / max(len(one_polys), len(other_polys)),
-    }
-
-
-@public
-def formula_similarity_fuzz(
-    one: Formula, other: Formula, curve: EllipticCurve, samples: int = 1000
-) -> Dict[str, float]:
-    if one.coordinate_model != other.coordinate_model:
-        warnings.warn("Mismatched coordinate model.")
-
-    output_matches = 0.0
-    iv_matches = 0.0
-    for _ in range(samples):
-        Paff = curve.affine_random()
-        Qaff = curve.affine_random()
-        Raff = curve.affine_add(Paff, Qaff)
-        P = Paff.to_model(one.coordinate_model, curve)
-        Q = Qaff.to_model(one.coordinate_model, curve)
-        R = Raff.to_model(one.coordinate_model, curve)
-        inputs = (P, Q, R)[: one.num_inputs]
-        with local(DefaultContext()) as ctx:
-            res_one = one(curve.prime, *inputs, **curve.parameters)
-        action_one = ctx.actions.get_by_index([0])
-        ivs_one = set(
-            map(attrgetter("value"), sum(action_one[0].intermediates.values(), []))
-        )
-        with local(DefaultContext()) as ctx:
-            res_other = other(curve.prime, *inputs, **curve.parameters)
-        action_other = ctx.actions.get_by_index([0])
-        ivs_other = set(
-            map(attrgetter("value"), sum(action_other[0].intermediates.values(), []))
-        )
-        iv_matches += len(ivs_one.intersection(ivs_other)) / max(
-            len(ivs_one), len(ivs_other)
-        )
-        one_coords = set(res_one)
-        other_coords = set(res_other)
-        output_matches += len(one_coords.intersection(other_coords)) / max(
-            len(one_coords), len(other_coords)
-        )
-    return {"output": output_matches / samples, "ivs": iv_matches / samples}
diff --git a/test/ec/test_formula.py b/test/ec/test_formula.py
index bb3d123..0e19d91 100644
--- a/test/ec/test_formula.py
+++ b/test/ec/test_formula.py
@@ -1,13 +1,33 @@
 import pickle
+from typing import Tuple
 
 import pytest
 from sympy import FF, symbols
-
+from importlib_resources import files, as_file
+import test.data.formulas
+from pyecsca.ec.formula.fliparoo import generate_fliparood_formulas
+from pyecsca.ec.formula.graph import rename_ivs
+from pyecsca.ec.formula.partitions import (
+    reduce_all_adds,
+    expand_all_muls,
+    expand_all_nopower2_muls,
+)
+from pyecsca.ec.formula.switch_sign import generate_switched_formulas
 from pyecsca.ec.mod import SymbolicMod, Mod
 from pyecsca.misc.cfg import TemporaryConfig
 from pyecsca.ec.error import UnsatisfiedAssumptionError
-from pyecsca.ec.params import get_params
+from pyecsca.ec.params import get_params, DomainParameters
 from pyecsca.ec.point import Point
+from pyecsca.ec.model import ShortWeierstrassModel, MontgomeryModel, TwistedEdwardsModel
+from pyecsca.ec.formula import (
+    AdditionEFDFormula,
+    DoublingEFDFormula,
+    LadderEFDFormula,
+    Formula,
+    AdditionFormula,
+    DoublingFormula,
+    LadderFormula,
+)
 
 
 @pytest.fixture()
@@ -62,39 +82,27 @@ def test_num_ops(add):
 
 
 def test_assumptions(secp128r1, mdbl):
-    res = mdbl(
-        secp128r1.curve.prime,
-        secp128r1.generator,
-        **secp128r1.curve.parameters
-    )
+    res = mdbl(secp128r1.curve.prime, secp128r1.generator, **secp128r1.curve.parameters)
     assert res is not None
 
-    coords = {
-        name: value * 5 for name, value in secp128r1.generator.coords.items()
-    }
+    coords = {name: value * 5 for name, value in secp128r1.generator.coords.items()}
     other = Point(secp128r1.generator.coordinate_model, **coords)
     with pytest.raises(UnsatisfiedAssumptionError):
-        mdbl(
-            secp128r1.curve.prime, other, **secp128r1.curve.parameters
-        )
+        mdbl(secp128r1.curve.prime, other, **secp128r1.curve.parameters)
     with TemporaryConfig() as cfg:
         cfg.ec.unsatisfied_formula_assumption_action = "ignore"
-        pt = mdbl(
-            secp128r1.curve.prime, other, **secp128r1.curve.parameters
-        )
+        pt = mdbl(secp128r1.curve.prime, other, **secp128r1.curve.parameters)
         assert pt is not None
 
 
 def test_parameters():
     jac_secp128r1 = get_params("secg", "secp128r1", "jacobian")
-    jac_dbl = jac_secp128r1.curve.coordinate_model.formulas[
-        "dbl-1998-hnm"
-    ]
+    jac_dbl = jac_secp128r1.curve.coordinate_model.formulas["dbl-1998-hnm"]
 
     res = jac_dbl(
         jac_secp128r1.curve.prime,
         jac_secp128r1.generator,
-        **jac_secp128r1.curve.parameters
+        **jac_secp128r1.curve.parameters,
     )
     assert res is not None
 
@@ -111,9 +119,7 @@ def test_symbolic(secp128r1, dbl):
         coords, **{key: SymbolicMod(symbols(key), p) for key in coords.variables}
     )
     symbolic_double = dbl(p, symbolic_point, **sympy_params)[0]
-    generator_double = dbl(
-        p, secp128r1.generator, **secp128r1.curve.parameters
-    )[0]
+    generator_double = dbl(p, secp128r1.generator, **secp128r1.curve.parameters)[0]
     for outer_var in coords.variables:
         symbolic_val = getattr(symbolic_double, outer_var).x
         generator_val = getattr(generator_double, outer_var).x
@@ -126,3 +132,279 @@ def test_symbolic(secp128r1, dbl):
 
 def test_pickle(add, dbl):
     assert add == pickle.loads(pickle.dumps(add))
+
+
+LIBRARY_FORMULAS = [
+    [
+        "add-bc-r1rv76-jac",
+        ShortWeierstrassModel,
+        "jacobian",
+        ("secg", "secp128r1"),
+        AdditionEFDFormula,
+    ],
+    [
+        "add-bc-r1rv76-mod",
+        ShortWeierstrassModel,
+        "modified",
+        ("secg", "secp128r1"),
+        AdditionEFDFormula,
+    ],
+    [
+        "dbl-bc-r1rv76-jac",
+        ShortWeierstrassModel,
+        "jacobian",
+        ("secg", "secp128r1"),
+        DoublingEFDFormula,
+    ],
+    [
+        "dbl-bc-r1rv76-mod",
+        ShortWeierstrassModel,
+        "modified",
+        ("secg", "secp128r1"),
+        DoublingEFDFormula,
+    ],
+    [
+        "dbl-bc-r1rv76-x25519",
+        MontgomeryModel,
+        "xz",
+        ("other", "Curve25519"),
+        DoublingEFDFormula,
+    ],
+    [
+        "ladd-bc-r1rv76-x25519",
+        MontgomeryModel,
+        "xz",
+        ("other", "Curve25519"),
+        LadderEFDFormula,
+    ],
+    [
+        "dbl-boringssl-p224",
+        ShortWeierstrassModel,
+        "jacobian-3",
+        ("secg", "secp224r1"),
+        DoublingEFDFormula,
+    ],
+    [
+        "add-boringssl-p224",
+        ShortWeierstrassModel,
+        "jacobian-3",
+        ("secg", "secp224r1"),
+        AdditionEFDFormula,
+    ],
+    [
+        "add-libressl-v382",
+        ShortWeierstrassModel,
+        "jacobian",
+        ("secg", "secp128r1"),
+        AdditionEFDFormula,
+    ],
+    [
+        "dbl-libressl-v382",
+        ShortWeierstrassModel,
+        "jacobian",
+        ("secg", "secp128r1"),
+        DoublingEFDFormula,
+    ],
+    [
+        "dbl-secp256k1-v040",
+        ShortWeierstrassModel,
+        "jacobian",
+        ("secg", "secp256k1"),
+        DoublingEFDFormula,
+    ],
+    [
+        "add-openssl-z256",
+        ShortWeierstrassModel,
+        "jacobian-3",
+        ("secg", "secp256r1"),
+        AdditionEFDFormula,
+    ],
+    [
+        "add-openssl-z256a",
+        ShortWeierstrassModel,
+        "jacobian-3",
+        ("secg", "secp256r1"),
+        AdditionEFDFormula,
+    ],
+    [
+        "ladd-openssl-x25519",
+        MontgomeryModel,
+        "xz",
+        ("other", "Curve25519"),
+        LadderEFDFormula,
+    ],
+    [
+        "ladd-hacl-x25519",
+        MontgomeryModel,
+        "xz",
+        ("other", "Curve25519"),
+        LadderEFDFormula,
+    ],
+    [
+        "dbl-hacl-x25519",
+        MontgomeryModel,
+        "xz",
+        ("other", "Curve25519"),
+        DoublingEFDFormula,
+    ],
+    [
+        "dbl-sunec-v21",
+        ShortWeierstrassModel,
+        "projective-3",
+        ("secg", "secp256r1"),
+        DoublingEFDFormula,
+    ],
+    [
+        "add-sunec-v21",
+        ShortWeierstrassModel,
+        "projective-3",
+        ("secg", "secp256r1"),
+        AdditionEFDFormula,
+    ],
+    [
+        "add-sunec-v21-ed25519",
+        TwistedEdwardsModel,
+        "extended",
+        ("other", "Ed25519"),
+        AdditionEFDFormula,
+    ],
+    [
+        "dbl-sunec-v21-ed25519",
+        TwistedEdwardsModel,
+        "extended",
+        ("other", "Ed25519"),
+        DoublingEFDFormula,
+    ],
+    [
+        "ladd-rfc7748",
+        MontgomeryModel,
+        "xz",
+        ("other", "Curve25519"),
+        LadderEFDFormula,
+    ],
+]
+
+
+@pytest.fixture(params=LIBRARY_FORMULAS)
+def library_formula_params(request) -> Tuple[Formula, DomainParameters]:
+    name, model, coords_name, param_spec, formula_type = request.param
+    model = model()
+    coordinate_model = model.coordinates[coords_name]
+    with as_file(files(test.data.formulas).joinpath(name)) as meta_path, as_file(
+        files(test.data.formulas).joinpath(name + ".op3")
+    ) as op3_path:
+        formula = formula_type(meta_path, op3_path, name, coordinate_model)
+    params = get_params(*param_spec, coords_name)
+    return formula, params
+
+
+def choose_curve(coordinate_model, name, library):
+    if library:
+        return next(filter(lambda x: x[0] == name, LIBRARY_FORMULAS))[3]
+    model = coordinate_model.curve_model
+    if model.__class__ == ShortWeierstrassModel:
+        return "secg", "secp128r1"
+    if model.__class__ == MontgomeryModel:
+        return "other", "Curve25519"
+    if model.__class__ == TwistedEdwardsModel:
+        return "other", "Ed25519"
+    raise NotImplementedError(model)
+
+
+def test_formula_graph(library_formula_params):
+    formula, params = library_formula_params
+    do_test_formula(rename_ivs(formula), params)
+
+
+def test_switch_sign(library_formula_params):
+    formula, params = library_formula_params
+    for switch_formula in generate_switched_formulas(formula):
+        do_test_formula(switch_formula, params)
+
+
+def test_fliparood_formula(library_formula_params):
+    formula, params = library_formula_params
+    for fliparood in generate_fliparood_formulas(formula):
+        do_test_formula(fliparood, params)
+
+
+def test_reductions(library_formula_params):
+    formula, params = library_formula_params
+    do_test_formula(reduce_all_adds(formula), params)
+
+
+def test_expansions(library_formula_params):
+    formula, params = library_formula_params
+    do_test_formula(expand_all_muls(formula), params)
+    do_test_formula(expand_all_nopower2_muls(formula), params)
+
+
+def do_test_formula(formula, params):
+    try:
+        do_test_formula0(formula, params)
+    except AssertionError:
+        print(formula.name)
+
+
+def do_test_formula0(formula, params):
+    coordinate_model = formula.coordinate_model
+    scale = coordinate_model.formulas.get("z", None)
+    if scale is None:
+        scale = coordinate_model.formulas.get("scale", None)
+
+    formula_type = formula.__class__
+    for _ in range(10):
+        Paff = params.curve.affine_random()
+        P2aff = params.curve.affine_double(Paff)
+        Qaff = params.curve.affine_random()
+        Q2aff = params.curve.affine_double(Qaff)
+        Raff = params.curve.affine_add(Paff, Qaff)
+        R2aff = params.curve.affine_double(Raff)
+        QRaff = params.curve.affine_add(Qaff, Raff)
+        P = Paff.to_model(coordinate_model, params.curve)
+        P2 = P2aff.to_model(coordinate_model, params.curve)
+        Q = Qaff.to_model(coordinate_model, params.curve)
+        Q2 = Q2aff.to_model(coordinate_model, params.curve)
+        R = Raff.to_model(coordinate_model, params.curve)
+        R2 = R2aff.to_model(coordinate_model, params.curve)  # noqa
+        QR = QRaff.to_model(coordinate_model, params.curve)
+        inputs = (P, Q, R)[: formula.num_inputs]
+        res = formula(params.curve.prime, *inputs, **params.curve.parameters)
+        if issubclass(formula_type, AdditionFormula):
+            try:
+                assert res[0].to_affine() == Raff
+            except NotImplementedError:
+                assert (
+                    scale(params.curve.prime, res[0], **params.curve.parameters)[0] == R
+                )
+        elif issubclass(formula_type, DoublingFormula):
+            try:
+                assert res[0].to_affine() == P2aff
+            except NotImplementedError:
+                assert (
+                    scale(params.curve.prime, res[0], **params.curve.parameters)[0]
+                    == P2
+                )
+        elif issubclass(formula_type, LadderFormula):
+            try:
+                assert res[0].to_affine() == Q2aff
+                assert res[1].to_affine() == QRaff
+            except NotImplementedError:
+                # print(scale(params.curve.prime, res[0], **params.curve.parameters)[0])
+                # print(scale(params.curve.prime, res[1], **params.curve.parameters)[0])
+                # print(P)
+                # print(Q)
+                # print(R)
+                # print(P2)
+                # print(Q2)
+                # print(R2)
+                # print(QR)
+                # print("------------------------------------")
+                assert (
+                    scale(params.curve.prime, res[1], **params.curve.parameters)[0]
+                    == QR
+                )
+                assert (
+                    scale(params.curve.prime, res[0], **params.curve.parameters)[0]
+                    == Q2
+                )
diff --git a/test/sca/test_structural.py b/test/sca/test_structural.py
index 970e4fc..2a883cc 100644
--- a/test/sca/test_structural.py
+++ b/test/sca/test_structural.py
@@ -11,7 +11,7 @@ from pyecsca.ec.formula import (
 )
 from pyecsca.ec.model import ShortWeierstrassModel, MontgomeryModel, TwistedEdwardsModel
 from pyecsca.ec.params import get_params
-from pyecsca.sca.re.structural import formula_similarity
+from pyecsca.ec.formula.metrics import formula_similarity
 
 
 def test_formula_similarity(secp128r1):