Merge branch 'J08nY:master' into CPA_correlations_tracking

25 files changed, 2593 insertions, 721 deletions

diff --git a/Makefile b/Makefile
index 537e885..5d0eaad 100644
--- a/Makefile
+++ b/Makefile
@@ -1,14 +1,3 @@
-EC_TESTS = test.ec.test_context test.ec.test_configuration test.ec.test_curve test.ec.test_formula \
-test.ec.test_params test.ec.test_key_agreement test.ec.test_key_generation test.ec.test_mod test.ec.test_model \
-test.ec.test_mult test.ec.test_naf test.ec.test_op test.ec.test_point test.ec.test_signature test.ec.test_transformations test.ec.test_regress \
-test.ec.test_divpoly
-
-SCA_TESTS = test.sca.test_align test.sca.test_combine test.sca.test_edit test.sca.test_filter test.sca.test_match test.sca.test_process \
-test.sca.test_sampling test.sca.test_target test.sca.test_test test.sca.test_trace test.sca.test_traceset test.sca.test_plot test.sca.test_rpa \
-test.sca.test_zvp test.sca.test_stacked_combine test.sca.test_leakage_models
-
-TESTS = ${EC_TESTS} ${SCA_TESTS}
-
 PERF_SCRIPTS = test.ec.perf_mod test.ec.perf_formula test.ec.perf_mult test.sca.perf_combine test.sca.perf_zvp
 
 test:
@@ -49,7 +38,7 @@ perf-plots:
 	python test/plots/plot_perf.py -d .perf
 
 doc-coverage:
-	interrogate -vv -nmps -e pyecsca/ec/std/.github/ -f 55 pyecsca
+	interrogate -c pyproject.toml -vv -nmps -f 55 pyecsca
 
 docs:
 	$(MAKE) -C docs apidoc
diff --git a/docs/references.rst b/docs/references.rst
index de3beae..675a13a 100644
--- a/docs/references.rst
+++ b/docs/references.rst
@@ -16,3 +16,5 @@ References
 .. [CO2002] Jean-Sébastien Coron. Resistance against Differential Power Analysis for Elliptic Curve Cryptosystems, https://link.springer.com/chapter/10.1007/3-540-48059-5_25
 .. [DJB02] D.J. Bernstein: Pippenger's Exponentiation Algorithm, https://cr.yp.to/papers/pippenger.pdf
 .. [SG14] Shay Gueron & Vlad Krasnov. Fast prime field elliptic-curve cryptography with 256-bit primes, https://link.springer.com/article/10.1007/s13389-014-0090-x
+.. [B51] Andrew D. Booth. A signed binary multiplication technique.
+.. [M61] O. L. Macsorley. High-speed arithmetic in binary computers.
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..3265131
--- /dev/null
+++ b/pyecsca/ec/formula/expand.py
@@ -0,0 +1,51 @@
+from typing import List, Callable, Any
+from public import public
+
+from . import Formula
+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
+
+
+@public
+def expand_formula_list(
+    formulas: List[EFDFormula], norm: Callable[[Formula], Any] = ivs_norm
+) -> List[EFDFormula]:
+    extended = reduce_with_similarity(formulas, norm)
+
+    fliparood: List[EFDFormula] = sum(list(map(recursive_fliparoo, extended)), [])
+    extended.extend(fliparood)
+    extended = reduce_with_similarity(extended, norm)
+
+    switch_signs: List[EFDFormula] = sum(
+        [list(generate_switched_formulas(f)) for f in extended], []
+    )
+    extended.extend(switch_signs)
+    extended = reduce_with_similarity(extended, norm)
+
+    extended.extend(list(map(reduce_all_adds, extended)))
+    extended = reduce_with_similarity(extended, norm)
+
+    extended.extend(list(map(expand_all_muls, extended)))
+    extended = reduce_with_similarity(extended, norm)
+
+    extended.extend(list(map(expand_all_nopower2_muls, extended)))
+    extended = reduce_with_similarity(extended, norm)
+
+    return extended
diff --git a/pyecsca/ec/formula/fliparoo.py b/pyecsca/ec/formula/fliparoo.py
new file mode 100644
index 0000000..c8d77ac
--- /dev/null
+++ b/pyecsca/ec/formula/fliparoo.py
@@ -0,0 +1,378 @@
+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:
+    """
+    Fliparoo is a chain of nodes N1->N2->...->Nk in EFDFormulaGraph for k>=2 such that:
+        - All Ni are * or All Ni are +/-
+        - For every Ni, except for Nk, the only outgoing node is Ni+1
+        - Neither of N1,...,Nk-1 is an output node
+    """
+
+    nodes: List[CodeOpNode]
+    graph: EFDFormulaGraph
+    operator: Optional[OpType]
+
+    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[CodeOpNode]):
+        for i, node in enumerate(nodes[:-1]):
+            if node.outgoing_nodes != [nodes[i + 1]]:
+                raise BadFliparoo
+            if node.output_node:
+                raise BadFliparoo
+
+    @property
+    def first(self):
+        return self.nodes[0]
+
+    @property
+    def last(self):
+        return self.nodes[-1]
+
+    def __len__(self):
+        return len(self.nodes)
+
+    def __repr__(self):
+        return "->".join(map(lambda x: x.__repr__(), self.nodes))
+
+    def previous(self, node: CodeOpNode) -> Optional[CodeOpNode]:
+        i = self.nodes.index(node)
+        if i == 0:
+            return None
+        return self.nodes[i - 1]
+
+    def __getitem__(self, i: int):
+        return self.nodes[i]
+
+    def __iter__(self):
+        return iter(self.nodes)
+
+    def __eq__(self, other):
+        return self.graph == other.graph and self.nodes == other.nodes
+
+    def deepcopy(self):
+        ngraph = self.graph.deepcopy()
+        nchain = [mirror_node(node, self.graph, ngraph) for node in self.nodes]
+        return self.__class__(nchain, ngraph)
+
+    def input_nodes(self) -> List[Node]:
+        input_nodes: List[Node] = []
+        for node in self:
+            input_nodes.extend(
+                filter(lambda x: x not in self.nodes, node.incoming_nodes)
+            )
+        return input_nodes
+
+    def slice(self, i: int, j: int):
+        return self.__class__(self.nodes[i:j], self.graph)
+
+
+class MulFliparoo(Fliparoo):
+    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:
+            raise BadFliparoo
+        self.operator = OpType.Mult
+
+
+class AddSubFliparoo(Fliparoo):
+    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]):
+            raise BadFliparoo
+
+
+class AddFliparoo(Fliparoo):
+    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:
+            raise BadFliparoo
+        self.operator = OpType.Add
+
+
+class BadFliparoo(Exception):
+    pass
+
+
+def find_fliparoos(
+    graph: EFDFormulaGraph, fliparoo_type: Optional[Type[Fliparoo]] = None
+) -> List[Fliparoo]:
+    """Finds a list of Fliparoos in a graph"""
+    fliparoos: List[Fliparoo] = []
+    for ilong_path in graph.find_all_paths():
+        long_path = ilong_path[1:]  # get rid of the input variables
+        fliparoo = largest_fliparoo(long_path, graph, fliparoo_type)  # type: ignore
+        if fliparoo and fliparoo not in fliparoos:
+            fliparoos.append(fliparoo)
+
+    # remove duplicities and fliparoos in inclusion
+    fliparoos = sorted(fliparoos, key=len, reverse=True)
+    longest_fliparoos: List[Fliparoo] = []
+    for fliparoo in fliparoos:
+        if not is_subfliparoo(fliparoo, longest_fliparoos):
+            longest_fliparoos.append(fliparoo)
+    return longest_fliparoos
+
+
+def is_subfliparoo(fliparoo: Fliparoo, longest_fliparoos: List[Fliparoo]) -> bool:
+    """Returns true if fliparoo is a part of any fliparoo in longest_fliparoos"""
+    for other_fliparoo in longest_fliparoos:
+        l1, l2 = len(fliparoo), len(other_fliparoo)
+        for i in range(l2 - l1):
+            if other_fliparoo.slice(i, i + l1) == fliparoo:
+                return True
+    return False
+
+
+def largest_fliparoo(
+    chain: List[CodeOpNode],
+    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):
+            subchain = chain[i : j + 1]
+            if fliparoo_type:
+                try:
+                    fliparoo_type(subchain, graph)
+                except BadFliparoo:
+                    continue
+            try:
+                return MulFliparoo(subchain, graph)
+            except BadFliparoo:
+                pass
+            try:
+                return AddSubFliparoo(subchain, graph)
+            except BadFliparoo:
+                pass
+    return None
+
+
+class SignedNode:
+
+    """
+    Represents a summand in an expression X1-X2+X3+X4-X5...
+    Used for creating +/- Fliparoos
+    """
+
+    node: CodeOpNode
+    sign: int
+
+    def __init__(self, node: CodeOpNode):
+        self.node = node
+        self.sign = 1
+
+    def __repr__(self):
+        s = {1: "+", -1: "-"}[self.sign]
+        return f"{s}{self.node.__repr__()}"
+
+
+class SignedSubGraph:
+    """Subgraph of an EFDFormula graph with signed nodes"""
+
+    def __init__(self, nodes: List[SignedNode], graph: EFDFormulaGraph):
+        self.nodes = nodes
+        self.supergraph = graph
+
+    def add_node(self, node: SignedNode):
+        self.nodes.append(node)
+
+    def remove_node(self, node: SignedNode):
+        self.nodes.remove(node)
+
+    def change_signs(self):
+        for node in self.nodes:
+            node.sign *= -1
+
+    def __getitem__(self, i):
+        return self.nodes[i]
+
+    @property
+    def components(self):
+        return len(self.nodes)
+
+    def deepcopy(self):
+        sgraph = self.supergraph.deepcopy()
+        return SignedSubGraph(
+            [mirror_node(n, self.supergraph, sgraph) for n in self.nodes], sgraph
+        )
+
+
+def mirror_node(node, graph, graph_copy):
+    """Finds the corresponding node in a copy of the graph"""
+    if isinstance(node, SignedNode):
+        ns = SignedNode(graph_copy.nodes[graph.node_index(node.node)])
+        ns.sign = node.sign
+        return ns
+    if isinstance(node, Node):
+        return graph_copy.nodes[graph.node_index(node)]
+    raise NotImplementedError
+
+
+class DummyNode(Node):
+    def __repr__(self):
+        return "Dummy node"
+
+    def label(self):
+        pass
+
+    def result(self):
+        pass
+
+
+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):
+            yield flip_graph.to_EFDFormula()
+
+
+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
+    signed_subgraphs = [signed_subgraph]
+    for _ in range(signed_subgraph.components - 1):
+        subgraphs_updated = []
+        for signed_subgraph in signed_subgraphs:
+            extended_subgraphs = combine_all_pairs_signed_nodes(
+                signed_subgraph, fliparoo
+            )
+            subgraphs_updated.extend(extended_subgraphs)
+        signed_subgraphs = subgraphs_updated
+
+    for signed_subgraph in signed_subgraphs:
+        graph = signed_subgraph.supergraph
+        assert signed_subgraph.components == 1
+        final_signed_node = signed_subgraph.nodes[0]
+        if final_signed_node.sign != 1:
+            continue
+        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))
+        reconnect_fliparoo_outputs(graph, final_node)
+        graph.update()
+        yield graph
+
+
+def extract_fliparoo_signed_inputs(fliparoo: Fliparoo) -> SignedSubGraph:
+    graph = fliparoo.graph
+    signed_inputs = SignedSubGraph([], graph)
+    for node in fliparoo:
+        prev = fliparoo.previous(node)
+        left, right = map(SignedNode, node.incoming_nodes)
+        if right.node != prev:
+            right.sign = -1 if node.is_sub else 1
+            signed_inputs.add_node(right)
+        if left.node != prev:
+            if node.is_sub and right.node == prev:
+                signed_inputs.change_signs()
+            signed_inputs.add_node(left)
+        if prev:
+            graph.remove_node(prev)
+    graph.remove_node(fliparoo.last)
+    return signed_inputs
+
+
+def disconnect_fliparoo_outputs(fliparoo: Fliparoo):
+    # remember positions of the last node with a DummyNode placeholder
+    dummy = DummyNode()
+    fliparoo.graph.add_node(dummy)
+    fliparoo.last.reconnect_outgoing_nodes(dummy)
+
+
+def reconnect_fliparoo_outputs(graph: EFDFormulaGraph, last_node: Node):
+    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))
+
+
+def combine_all_pairs_signed_nodes(
+    signed_subgraph: SignedSubGraph, fliparoo: Fliparoo
+) -> List[SignedSubGraph]:
+    signed_subgraphs = []
+    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)
+            signed_subgraphs.append(csigned_subgraph)
+    return signed_subgraphs
+
+
+def combine_signed_nodes(
+    subgraph: SignedSubGraph,
+    left_signed_node: SignedNode,
+    right_signed_node: SignedNode,
+    fliparoo: Fliparoo,
+):
+    left_node, right_node = left_signed_node.node, right_signed_node.node
+    sign = 1
+    operator = OpType.Mult
+    if isinstance(fliparoo, AddSubFliparoo):
+        s0, s1 = left_signed_node.sign, right_signed_node.sign
+        if s0 == 1:
+            operator = OpType.Add if s1 == 1 else OpType.Sub
+
+        if s0 == -1 and s1 == 1:
+            operator = OpType.Sub
+            left_node, right_node = right_node, left_node
+
+        # propagate the sign
+        if s0 == -1 and s1 == -1:
+            operator = OpType.Add
+            sign = -1
+
+    new_node = CodeOpNode.from_str(
+        f"Fliparoo{id(left_node)}_{id(operator)}_{id(sign)}_{id(right_node)}", left_node.result, operator, right_node.result
+    )
+    new_node.incoming_nodes = [left_node, right_node]
+    left_node.outgoing_nodes.append(new_node)
+    right_node.outgoing_nodes.append(new_node)
+    subgraph.supergraph.add_node(new_node)
+    new_node = SignedNode(new_node)
+    new_node.sign = sign
+    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/graph.py b/pyecsca/ec/formula/graph.py
new file mode 100644
index 0000000..1473858
--- /dev/null
+++ b/pyecsca/ec/formula/graph.py
@@ -0,0 +1,436 @@
+from .efd import (
+    EFDFormula,
+    DoublingEFDFormula,
+    AdditionEFDFormula,
+    LadderEFDFormula,
+    DifferentialAdditionEFDFormula,
+)
+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, Optional, MutableMapping, Any
+from copy import deepcopy
+from abc import ABC, abstractmethod
+
+
+class Node(ABC):
+    def __init__(self):
+        self.incoming_nodes = []
+        self.outgoing_nodes = []
+        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
+
+    def reconnect_outgoing_nodes(self, destination):
+        destination.output_node = self.output_node
+        for out in self.outgoing_nodes:
+            out.incoming_nodes = [
+                n if n != self else destination for n in out.incoming_nodes
+            ]
+            destination.outgoing_nodes.append(out)
+
+
+class ConstantNode(Node):
+    color = "#b41f44"
+
+    def __init__(self, i: int):
+        super().__init__()
+        self.value = i
+
+    @property
+    def label(self) -> str:
+        return str(self.value)
+
+    @property
+    def result(self) -> str:
+        return str(self.value)
+
+    def __repr__(self) -> str:
+        return f"Node({self.value})"
+
+
+class CodeOpNode(Node):
+    color = "#1f78b4"
+
+    def __init__(self, op: CodeOp):
+        super().__init__()
+        self.op = op
+        assert self.op.operator in [
+            OpType.Sub,
+            OpType.Add,
+            OpType.Id,
+            OpType.Sqr,
+            OpType.Mult,
+            OpType.Div,
+            OpType.Pow,
+            OpType.Inv,
+            OpType.Neg,
+        ], self.op.operator
+
+    @classmethod
+    def from_str(cls, result: str, left, operator, right):
+        opstr = f"{result} = {left if left is not None else ''}{operator.op_str}{right if right is not None else ''}"
+        return cls(CodeOp(parse(opstr.replace("^", "**"))))
+
+    @property
+    def label(self) -> str:
+        return f"{self.op.result}:{self.op.operator.op_str}"
+
+    @property
+    def result(self) -> str:
+        return str(self.op.result)
+
+    @property
+    def optype(self) -> OpType:
+        return self.op.operator
+
+    @property
+    def is_sub(self) -> bool:
+        return self.optype == OpType.Sub
+
+    @property
+    def is_mul(self) -> bool:
+        return self.optype == OpType.Mult
+
+    @property
+    def is_add(self) -> bool:
+        return self.optype == OpType.Add
+
+    @property
+    def is_id(self) -> bool:
+        return self.optype == OpType.Id
+
+    @property
+    def is_sqr(self) -> bool:
+        return self.optype == OpType.Sqr
+
+    @property
+    def is_pow(self) -> bool:
+        return self.optype == OpType.Pow
+
+    @property
+    def is_inv(self) -> bool:
+        return self.optype == OpType.Inv
+
+    @property
+    def is_div(self) -> bool:
+        return self.optype == OpType.Div
+
+    @property
+    def is_neg(self) -> bool:
+        return self.optype == OpType.Neg
+
+    def __repr__(self) -> str:
+        return f"Node({self.op})"
+
+
+class InputNode(Node):
+    color = "#b41f44"
+
+    def __init__(self, input: str):
+        super().__init__()
+        self.input = input
+        self.input_node = True
+
+    @property
+    def label(self) -> str:
+        return self.input
+
+    @property
+    def result(self) -> str:
+        return self.input
+
+    def __repr__(self) -> str:
+        return f"Node({self.input})"
+
+
+def formula_input_variables(formula: EFDFormula) -> List[str]:
+    return (
+        list(formula.inputs)
+        + formula.parameters
+        + formula.coordinate_model.curve_model.parameter_names
+    )
+
+
+# temporary solution
+class ModifiedEFDFormula(EFDFormula):
+    def __eq__(self, other):
+        if not isinstance(other, ModifiedEFDFormula):
+            return False
+        return (
+            self.name == other.name and self.coordinate_model == other.coordinate_model and self.code == other.code
+        )
+
+
+class ModifiedDoublingEFDFormula(DoublingEFDFormula, ModifiedEFDFormula):
+    pass
+
+
+class ModifiedAdditionEFDFormula(AdditionEFDFormula, ModifiedEFDFormula):
+    pass
+
+
+class ModifiedDifferentialAdditionEFDFormula(
+    DifferentialAdditionEFDFormula, ModifiedEFDFormula
+):
+    pass
+
+
+class ModifiedLadderEFDFormula(LadderEFDFormula, ModifiedEFDFormula):
+    pass
+
+
+class EFDFormulaGraph:
+    nodes: List[Node]
+    input_nodes: MutableMapping[str, InputNode]
+    output_names: Set[str]
+    roots: List[Node]
+    coordinate_model: Any
+
+    def __init__(self, formula: EFDFormula, rename=True):
+        self._formula = formula  # TODO remove, its here only for to_EFDFormula
+        self.coordinate_model = formula.coordinate_model
+        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: Dict[str, Node] = self.input_nodes.copy()  # type: ignore
+        constants: Dict[int, Node] = {}
+        for op in formula.code:
+            code_node = CodeOpNode(op)
+            for side in (op.left, op.right):
+                if side is None:
+                    continue
+                if isinstance(side, int):
+                    if side in constants:
+                        parent_node = constants[side]
+                    else:
+                        parent_node = ConstantNode(side)
+                    self.nodes.append(parent_node)
+                    self.roots.append(parent_node)
+                else:
+                    parent_node = discovered_nodes[side]
+                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 node not in self.roots:
+                self.roots.append(node)
+        if rename:
+            self.reindex()
+
+    def node_index(self, node: Node) -> int:
+        return self.nodes.index(node)
+
+    def deepcopy(self):
+        return deepcopy(self)
+
+    def to_EFDFormula(self) -> ModifiedEFDFormula:
+        # TODO rewrite
+        new_graph = deepcopy(self)
+        new_formula = new_graph._formula
+        new_formula.code = list(
+            map(
+                lambda x: x.op,  # type: ignore
+                filter(lambda n: n not in new_graph.roots, new_graph.nodes),
+            )
+        )
+        casting = {
+            AdditionEFDFormula: ModifiedAdditionEFDFormula,
+            DoublingEFDFormula: ModifiedDoublingEFDFormula,
+            DifferentialAdditionEFDFormula: ModifiedDifferentialAdditionEFDFormula,
+            LadderEFDFormula: ModifiedLadderEFDFormula,
+        }
+        if new_formula.__class__ not in set(casting.values()):
+            new_formula.__class__ = casting[new_formula.__class__]
+        return new_formula  # type: ignore
+
+    def networkx_graph(self) -> nx.DiGraph:
+        graph = nx.DiGraph()
+        vertices = list(range(len(self.nodes)))
+        graph.add_nodes_from(vertices)
+        stack = self.roots.copy()
+        while stack:
+            node = stack.pop()
+            for out in node.outgoing_nodes:
+                stack.append(out)
+                graph.add_edge(self.node_index(node), self.node_index(out))
+        return graph
+
+    def levels(self) -> List[List[Node]]:
+        stack = self.roots.copy()
+        levels = [(n, 0) for n in stack]
+        level_counter = 1
+        while stack:
+            tmp_stack = []
+            while stack:
+                node = stack.pop()
+                levels.append((node, level_counter))
+                for out in node.outgoing_nodes:
+                    tmp_stack.append(out)
+            stack = tmp_stack
+            level_counter += 1
+        # separate into lists
+
+        level_lists: List[List[Node]] = [[] for _ in range(level_counter)]
+        discovered = []
+        for node, l in reversed(levels):
+            if node not in discovered:
+                level_lists[l].append(node)
+                discovered.append(node)
+        return level_lists
+
+    def output_nodes(self) -> List[Node]:
+        return list(filter(lambda x: x.output_node, self.nodes))
+
+    def planar_positions(self) -> Dict[int, Tuple[float, float]]:
+        positions = {}
+        for i, level in enumerate(self.levels()):
+            for j, node in enumerate(level):
+                positions[self.node_index(node)] = (
+                    0.1 * float(i**2) + 3 * float(j),
+                    -6 * float(i),
+                )
+        return positions
+
+    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)
+        colors = [self.nodes[n].color for n in gnx.nodes]
+        labels = {n: self.nodes[n].label for n in gnx.nodes}
+        nx.draw(gnx, pos, node_color=colors, node_size=2000, labels=labels)
+        if filename:
+            plt.savefig(filename)
+            plt.close()
+        else:
+            plt.show()
+
+    def find_all_paths(self) -> List[List[Node]]:
+        gnx = self.networkx_graph()
+        index_paths = []
+        for u in self.roots:
+            iu = self.node_index(u)
+            for v in self.output_nodes():
+                iv = self.node_index(v)
+                index_paths.extend(nx.all_simple_paths(gnx, iu, iv))
+        paths = []
+        for p in index_paths:
+            paths.append([self.nodes[v] for v in p])
+        return paths
+
+    def reorder(self):
+        self.nodes = sum(self.levels(), [])
+
+    def remove_node(self, node):
+        self.nodes.remove(node)
+        if node in self.roots:
+            self.roots.remove(node)
+        for in_node in node.incoming_nodes:
+            in_node.outgoing_nodes = list(
+                filter(lambda x: x != node, in_node.outgoing_nodes)
+            )
+        for out_node in node.outgoing_nodes:
+            out_node.incoming_nodes = list(
+                filter(lambda x: x != node, out_node.incoming_nodes)
+            )
+
+    def add_node(self, node):
+        if isinstance(node, ConstantNode):
+            self.roots.append(node)
+        self.nodes.append(node)
+
+    def reindex(self):
+        results: Dict[str, str] = {}
+        counter = 0
+        for node in self.nodes:
+            if not isinstance(node, CodeOpNode):
+                continue
+            op = node.op
+            result, left, operator, right = (
+                op.result,
+                op.left,
+                op.operator.op_str,
+                op.right,
+            )
+            if left in results:
+                left = results[left]
+            if right in results:
+                right = results[right]
+            if not node.output_node:
+                new_result = f"iv{counter}"
+                counter += 1
+            else:
+                new_result = result
+            opstr = f"{new_result} = {left if left is not None else ''}{operator}{right if right is not None else ''}"
+            results[result] = new_result
+            node.op = CodeOp(parse(opstr.replace("^", "**")))
+
+    def update(self):
+        self.reorder()
+        self.reindex()
+
+    def print(self):
+        for node in self.nodes:
+            print(node)
+
+
+def rename_ivs(formula: EFDFormula):
+    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/partitions.py b/pyecsca/ec/formula/partitions.py
new file mode 100644
index 0000000..9ea108c
--- /dev/null
+++ b/pyecsca/ec/formula/partitions.py
@@ -0,0 +1,378 @@
+from typing import List, Any, Generator
+from ast import parse
+from ..op import OpType, CodeOp
+from .graph import (
+    EFDFormulaGraph,
+    CodeOpNode,
+    ConstantNode,
+    Node,
+)
+from .fliparoo import find_fliparoos, AddFliparoo, MulFliparoo
+from copy import deepcopy
+from .efd import EFDFormula
+
+
+def reduce_all_adds(formula: EFDFormula, rename=True) -> EFDFormula:
+    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)
+    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(formula, rename)
+    enodes = find_expansion_nodes(graph)
+    for enode in enodes:
+        expand_mul(graph, enode, copy=False)
+    return graph.to_EFDFormula()
+
+
+def expand_all_nopower2_muls(formula: EFDFormula, rename=True) -> EFDFormula:
+    graph = EFDFormulaGraph(formula, rename)
+    enodes = find_expansion_nodes(graph, nopower2=True)
+    for enode in enodes:
+        expand_mul(graph, enode, copy=False)
+    return graph.to_EFDFormula()
+
+
+def find_single_input_add_fliparoos(graph: EFDFormulaGraph) -> List[AddFliparoo]:
+    fliparoos = find_fliparoos(graph, AddFliparoo)
+    single_input_fliparoos = []
+    for fliparoo in fliparoos:
+        found = False
+        for i in range(len(fliparoo), 1, -1):
+            subfliparoo = fliparoo.slice(0, i)
+            if len(set(subfliparoo.input_nodes())) == 1:
+                found = True
+                break
+        if found:
+            s = subfliparoo.slice(0, i)
+            single_input_fliparoos.append(s)
+    return single_input_fliparoos
+
+
+def find_constant_mul_fliparoos(graph: EFDFormulaGraph) -> List[MulFliparoo]:
+    fliparoos = find_fliparoos(graph, MulFliparoo)
+    constant_mul_fliparoo = []
+    for fliparoo in fliparoos:
+        found = False
+        for i in range(len(fliparoo), 1, -1):
+            subfliparoo = fliparoo.slice(0, i)
+            nonconstant_inputs = list(
+                filter(
+                    lambda x: not isinstance(x, ConstantNode), subfliparoo.input_nodes()
+                )
+            )
+            if len(nonconstant_inputs) != 1:
+                continue
+            inode = nonconstant_inputs[0]
+            if inode not in fliparoo.first.incoming_nodes:
+                continue
+            if not sum(
+                1
+                for node in fliparoo.first.incoming_nodes
+                if isinstance(node, ConstantNode)
+            ):
+                continue
+            found = True
+            break
+        if found:
+            s = subfliparoo.slice(0, i)
+            constant_mul_fliparoo.append(s)
+    return constant_mul_fliparoo
+
+
+def find_expansion_nodes(graph: EFDFormulaGraph, nopower2=False) -> List[Node]:
+    expansion_nodes: List[Node] = []
+    for node in graph.nodes:
+        if not isinstance(node, CodeOpNode) or not node.is_mul:
+            continue
+        for par in node.incoming_nodes:
+            if isinstance(par, ConstantNode):
+                if nopower2 and is_power_of_2(par.value):
+                    continue
+                expansion_nodes.append(node)
+                break
+    return expansion_nodes
+
+
+def is_power_of_2(n: int) -> bool:
+    while n > 1:
+        if n & 1 == 1:
+            return False
+        n >>= 1
+    return True
+
+
+def reduce_all_XplusX(graph: EFDFormulaGraph):
+    adds = find_all_XplusX(graph)
+    for node in adds:
+        reduce_XplusX(graph, node)
+    graph.update()
+
+
+def find_all_XplusX(graph) -> List[CodeOpNode]:
+    adds = []
+    for node in graph.nodes:
+        if not isinstance(node, CodeOpNode) or not node.is_add:
+            continue
+        if node.incoming_nodes[0] == node.incoming_nodes[1]:
+            adds.append(node)
+    return adds
+
+
+def reduce_XplusX(graph: EFDFormulaGraph, node: CodeOpNode):
+    inode = node.incoming_nodes[0]
+    const_node = ConstantNode(2)
+    node.incoming_nodes[1] = const_node
+    const_node.outgoing_nodes = [node]
+    graph.add_node(const_node)
+    inode.outgoing_nodes = list(filter(lambda x: x != node, inode.outgoing_nodes))
+    inode.outgoing_nodes.append(node)
+    opstr = f"{node.result} = {inode.result}{OpType.Mult.op_str}{const_node.value}"
+    node.op = CodeOp(parse(opstr))
+
+
+def reduce_mul_fliparoo(fliparoo: MulFliparoo, copy=True):
+    if copy:
+        fliparoo = fliparoo.deepcopy()
+
+    first, last = fliparoo.first, fliparoo.last
+    inode = next(
+        filter(lambda x: not isinstance(x, ConstantNode), first.incoming_nodes)
+    )
+    const_nodes: List[ConstantNode] = [node for node in fliparoo.input_nodes() if isinstance(node, ConstantNode)]
+    sum_const_node = ConstantNode(sum(v.value for v in const_nodes))
+    fliparoo.graph.add_node(sum_const_node)
+
+    inode.outgoing_nodes = [n if n != first else last for n in inode.outgoing_nodes]
+    last.incoming_nodes = [inode, sum_const_node]
+    sum_const_node.outgoing_nodes = [last]
+
+    opstr = f"{last.result} = {inode.result}{OpType.Mult.op_str}{sum_const_node.value}"
+    last.op = CodeOp(parse(opstr))
+
+    for node in fliparoo:
+        if node == last:
+            continue
+        fliparoo.graph.remove_node(node)
+
+    for node in const_nodes:
+        if not node.outgoing_nodes:
+            fliparoo.graph.remove_node(node)
+
+    fliparoo.graph.update()
+
+    return fliparoo.graph
+
+
+def reduce_add_fliparoo(fliparoo: AddFliparoo, copy=True):
+    if copy:
+        fliparoo = fliparoo.deepcopy()
+    first, last = fliparoo.first, fliparoo.last
+    par = first.incoming_nodes[0]
+    const_node = ConstantNode(len(fliparoo) + 1)
+    fliparoo.graph.add_node(const_node)
+    mul_node = CodeOpNode.from_str(
+        last.result, const_node.result, OpType.Mult, par.result
+    )
+    fliparoo.graph.add_node(mul_node)
+    mul_node.incoming_nodes = [const_node, par]
+    par.outgoing_nodes.append(mul_node)
+    const_node.outgoing_nodes.append(mul_node)
+    mul_node.output_node = last.output_node
+    last.reconnect_outgoing_nodes(mul_node)
+    for node in fliparoo:
+        fliparoo.graph.remove_node(node)
+
+    fliparoo.graph.update()
+
+    return fliparoo.graph
+
+
+def expand_mul(graph: EFDFormulaGraph, node: Node, copy=True) -> EFDFormulaGraph:
+    if copy:
+        i = graph.node_index(node)
+        graph = deepcopy(graph)
+        node = graph.nodes[i]
+
+    const_par = next(filter(lambda x: isinstance(x, ConstantNode), node.incoming_nodes))
+    par = next(filter(lambda x: not isinstance(x, ConstantNode), node.incoming_nodes))
+    initial_node = CodeOpNode.from_str(node.result, par.result, OpType.Add, par.result)
+    graph.add_node(initial_node)
+    initial_node.incoming_nodes = [par, par]
+    par.outgoing_nodes.extend([initial_node, initial_node])
+    prev_node = initial_node
+    for _ in range(const_par.value - 2):
+        anode = CodeOpNode.from_str(
+            node.result, prev_node.result, OpType.Add, par.result
+        )
+        anode.incoming_nodes = [prev_node, par]
+        par.outgoing_nodes.append(anode)
+        graph.add_node(anode)
+        prev_node.outgoing_nodes = [anode]
+        prev_node = anode
+
+    prev_node.output_node = node.output_node
+    node.reconnect_outgoing_nodes(prev_node)
+    graph.remove_node(node)
+    graph.remove_node(const_par)
+    graph.update()
+
+    return graph
+
+
+class Partition:
+    value: int
+    parts: List["Partition"]
+
+    def __init__(self, n: int):
+        self.value = n
+        self.parts = []
+
+    @property
+    def is_final(self):
+        return not self.parts
+
+    def __repr__(self):
+        if self.is_final:
+            return f"({self.value})"
+        l, r = self.parts
+        return f"({l.__repr__()},{r.__repr__()})"
+
+    def __add__(self, other):
+        a = Partition(self.value + other.value)
+        a.parts = [self, other]
+        return a
+
+    def __eq__(self, other):
+        if self.value != other.value:
+            return False
+        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)
+
+
+def compute_partitions(n: int) -> List[Partition]:
+    partitions = [Partition(n)]
+    for d in range(1, n // 2 + 1):
+        n_d = n - d
+        for partition_dp in compute_partitions(d):
+            for partition_n_dp in compute_partitions(n_d):
+                partitions.append(partition_dp + partition_n_dp)
+    # remove duplicates
+    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(formula, rename)
+    enodes = find_expansion_nodes(graph)
+    for enode in enodes:
+        for part_graph in generate_all_node_partitions(graph, enode):
+            yield part_graph.to_EFDFormula()
+
+
+def generate_all_node_partitions(
+    original_graph: EFDFormulaGraph, node: Node
+) -> Generator[EFDFormulaGraph, Any, None]:
+    const_par = next(filter(lambda x: isinstance(x, ConstantNode), node.incoming_nodes))
+    const_par_value = const_par.value
+
+    par = next(filter(lambda x: not isinstance(x, ConstantNode), node.incoming_nodes))
+    i, ic, ip = (
+        original_graph.node_index(node),
+        original_graph.node_index(const_par),
+        original_graph.node_index(par),
+    )
+
+    for partition in compute_partitions(const_par_value):
+        if partition.is_final:
+            continue
+
+        # copy
+        graph = deepcopy(original_graph)
+        node, const_par, par = graph.nodes[i], graph.nodes[ic], graph.nodes[ip]
+        graph.remove_node(const_par)
+        lresult, rresult = f"{node.result}L", f"{node.result}R"
+        empty_left_node = CodeOpNode.from_str(lresult, "PART", OpType.Add, "PART")
+        empty_right_node = CodeOpNode.from_str(rresult, "PART", OpType.Add, "PART")
+        addition_node = CodeOpNode.from_str(node.result, lresult, OpType.Add, rresult)
+        graph.add_node(empty_left_node)
+        graph.add_node(empty_right_node)
+        graph.add_node(addition_node)
+        addition_node.outgoing_nodes = node.outgoing_nodes
+        addition_node.output_node = node.output_node
+        addition_node.incoming_nodes = [empty_left_node, empty_right_node]
+        empty_left_node.outgoing_nodes = [addition_node]
+        empty_right_node.outgoing_nodes = [addition_node]
+
+        left, right = partition.parts
+        partition_node(graph, empty_left_node, left, par)
+        partition_node(graph, empty_right_node, right, par)
+
+        graph.remove_node(node)
+        graph.update()
+        yield graph
+
+
+def partition_node(
+    graph: EFDFormulaGraph, node: CodeOpNode, partition: Partition, source_node: Node
+):
+    if partition.is_final and partition.value == 1:
+        # source node will take the role of node
+        # note: node has precisely one output node, since it was created during partitions
+        assert len(node.outgoing_nodes) == 1
+        child = node.outgoing_nodes[0]
+        source_node.outgoing_nodes.append(child)
+
+        left, right = child.incoming_nodes[0].result, child.incoming_nodes[1].result
+        if child.incoming_nodes[0] == node:
+            left = source_node.result
+            child.incoming_nodes[0] = source_node
+        else:
+            right = source_node.result
+            child.incoming_nodes[1] = source_node
+        opstr = f"{child.result} = {left}{child.optype.op_str}{right}"
+        child.op = CodeOp(parse(opstr))
+        graph.remove_node(node)
+        return
+
+    if partition.is_final:
+        source_node.outgoing_nodes.append(node)
+        const_node = ConstantNode(partition.value)
+        graph.add_node(const_node)
+        opstr = (
+            f"{node.result} = {source_node.result}{OpType.Mult.op_str}{partition.value}"
+        )
+        node.op = CodeOp(parse(opstr))
+        node.incoming_nodes = [source_node, const_node]
+        const_node.outgoing_nodes = [node]
+        return
+
+    lresult, rresult = f"{node.result}L", f"{node.result}R"
+    empty_left_node = CodeOpNode.from_str(lresult, "PART", OpType.Add, "PART")
+    empty_right_node = CodeOpNode.from_str(rresult, "PART", OpType.Add, "PART")
+
+    opstr = f"{node.result} = {lresult}{OpType.Add.op_str}{rresult}"
+    node.op = CodeOp(parse(opstr))
+    graph.add_node(empty_left_node)
+    graph.add_node(empty_right_node)
+
+    node.incoming_nodes = [empty_left_node, empty_right_node]
+    empty_left_node.outgoing_nodes = [node]
+    empty_right_node.outgoing_nodes = [node]
+
+    left, right = partition.parts
+    partition_node(graph, empty_left_node, left, source_node)
+    partition_node(graph, empty_right_node, right, source_node)
diff --git a/pyecsca/ec/formula/switch_sign.py b/pyecsca/ec/formula/switch_sign.py
new file mode 100644
index 0000000..1acef5b
--- /dev/null
+++ b/pyecsca/ec/formula/switch_sign.py
@@ -0,0 +1,131 @@
+from typing import Dict, Iterator, List, Any
+from ast import parse
+from ..op import OpType, CodeOp
+from .graph import EFDFormulaGraph, ConstantNode, Node, CodeOpNode
+from itertools import chain, combinations
+from .efd import EFDFormula
+from ..point import Point
+from ..mod import Mod
+
+
+def generate_switched_formulas(
+    formula: EFDFormula, rename=True
+) -> Iterator[EFDFormula]:
+    graph = EFDFormulaGraph(formula, rename)
+    for node_combination in subnode_lists(graph):
+        try:
+            yield switch_sign(graph, node_combination).to_EFDFormula()
+        except BadSignSwitch:
+            continue
+
+
+def subnode_lists(graph: EFDFormulaGraph):
+    return powerlist(filter(lambda x: x not in graph.roots and x.is_sub, graph.nodes))
+
+
+def switch_sign(graph: EFDFormulaGraph, node_combination) -> EFDFormulaGraph:
+    nodes_i = [graph.node_index(node) for node in node_combination]
+    graph = graph.deepcopy()
+    node_combination = set(graph.nodes[node_i] for node_i in nodes_i)
+    output_signs = {out: 1 for out in graph.output_names}
+
+    queue = []
+    for node in node_combination:
+        change_sides(node)
+        if node.output_node:
+            output_signs[node.result] = -1
+        queue.extend([(out, node.result) for out in node.outgoing_nodes])
+
+    while queue:
+        node, variable = queue.pop()
+        queue = switch_sign_propagate(node, variable, output_signs) + queue
+
+    sign_test(output_signs, graph.coordinate_model)
+    return graph
+
+
+def sign_test(output_signs: Dict[str, int], coordinate_model: Any):
+    scale = coordinate_model.formulas.get("z", None)
+    if scale is None:
+        scale = coordinate_model.formulas.get("scale", None)
+    p = 7
+    out_inds = set(map(lambda x: "".join([o for o in x if o.isdigit()]), output_signs))
+    for ind in out_inds:
+        point_dict = {}
+        for out, sign in output_signs.items():
+            if not out.endswith(ind):
+                continue
+            out_var = out[:out.index(ind)]
+            if not out_var.isalpha():
+                continue
+            point_dict[out_var] = Mod(sign, p)
+        point = Point(coordinate_model, **point_dict)
+        try:
+            apoint = point.to_affine()
+        except NotImplementedError:
+            apoint = scale(p, point)[0]
+        if set(apoint.coords.values()) != set([Mod(1, p)]):
+            raise BadSignSwitch
+
+
+class BadSignSwitch(Exception):
+    pass
+
+
+def switch_sign_propagate(
+    node: CodeOpNode, variable: str, output_signs: Dict[str, int]
+):
+    if node.is_add:
+        if variable == node.incoming_nodes[1].result:
+            node.op = change_operator(node.op, OpType.Sub)
+            return []
+        change_sides(node)
+        node.op = change_operator(node.op, OpType.Sub)
+        return []
+    if node.is_id or node.is_neg:
+        output_signs[node.result] *= -1
+        return [(child, node.result) for child in node.outgoing_nodes]
+
+    if node.is_sqr:
+        return []
+    if node.is_sub:
+        if node.incoming_nodes[0].result == variable:
+            node.op = change_operator(node.op, OpType.Add)
+            if node.output_node:
+                output_signs[node.result] *= -1
+            return [(child, node.result) for child in node.outgoing_nodes]
+        node.op = change_operator(node.op, OpType.Add)
+        return []
+    if node.is_pow:
+        exponent = next(
+            filter(lambda n: isinstance(n, ConstantNode), node.incoming_nodes)
+        )
+        if exponent.value % 2 == 0:
+            return []
+    if node.output_node:
+        output_signs[node.result] *= -1
+    assert node.is_mul or node.is_pow or node.is_inv or node.is_div
+    return [(child, node.result) for child in node.outgoing_nodes]
+
+
+def change_operator(op, new_operator):
+    result, left, right = op.result, op.left, op.right
+    opstr = f"{result} = {left if left is not None else ''}{new_operator.op_str}{right if right is not None else ''}"
+    return CodeOp(parse(opstr.replace("^", "**")))
+
+
+def change_sides(node):
+    op = node.op
+    result, left, operator, right = op.result, op.left, op.operator.op_str, op.right
+    left, right = right, left
+    opstr = f"{result} = {left if left is not None else ''}{operator}{right if right is not None else ''}"
+    node.op = CodeOp(parse(opstr.replace("^", "**")))
+    node.incoming_nodes[1], node.incoming_nodes[0] = (
+        node.incoming_nodes[0],
+        node.incoming_nodes[1],
+    )
+
+
+def powerlist(iterable: Iterator) -> List:
+    s = list(iterable)
+    return list(chain.from_iterable(combinations(s, r) for r in range(len(s) + 1)))
diff --git a/pyecsca/ec/mult/window.py b/pyecsca/ec/mult/window.py
index f85a58a..d025cc1 100644
--- a/pyecsca/ec/mult/window.py
+++ b/pyecsca/ec/mult/window.py
@@ -4,15 +4,22 @@ from typing import Optional, MutableMapping
 from public import public
 
 from ..params import DomainParameters
-from .base import ScalarMultiplier, AccumulationOrder, ScalarMultiplicationAction, PrecomputationAction, \
-    ProcessingDirection, AccumulatorMultiplier
+from .base import (
+    ScalarMultiplier,
+    AccumulationOrder,
+    ScalarMultiplicationAction,
+    PrecomputationAction,
+    ProcessingDirection,
+    AccumulatorMultiplier,
+)
 from ..formula import (
     AdditionFormula,
     DoublingFormula,
     ScalingFormula,
+    NegationFormula,
 )
 from ..point import Point
-from ..scalar import convert_base, sliding_window_rtl, sliding_window_ltr
+from ..scalar import convert_base, sliding_window_rtl, sliding_window_ltr, booth_window
 
 
 @public
@@ -34,17 +41,21 @@ class SlidingWindowMultiplier(AccumulatorMultiplier, ScalarMultiplier):
     _points: MutableMapping[int, Point]
 
     def __init__(
-            self,
-            add: AdditionFormula,
-            dbl: DoublingFormula,
-            width: int,
-            scl: Optional[ScalingFormula] = None,
-            recoding_direction: ProcessingDirection = ProcessingDirection.LTR,
-            accumulation_order: AccumulationOrder = AccumulationOrder.PeqPR,
-            short_circuit: bool = True,
+        self,
+        add: AdditionFormula,
+        dbl: DoublingFormula,
+        width: int,
+        scl: Optional[ScalingFormula] = None,
+        recoding_direction: ProcessingDirection = ProcessingDirection.LTR,
+        accumulation_order: AccumulationOrder = AccumulationOrder.PeqPR,
+        short_circuit: bool = True,
     ):
         super().__init__(
-            short_circuit=short_circuit, accumulation_order=accumulation_order, add=add, dbl=dbl, scl=scl
+            short_circuit=short_circuit,
+            accumulation_order=accumulation_order,
+            add=add,
+            dbl=dbl,
+            scl=scl,
         )
         self.width = width
         self.recoding_direction = recoding_direction
@@ -55,7 +66,13 @@ class SlidingWindowMultiplier(AccumulatorMultiplier, ScalarMultiplier):
     def __eq__(self, other):
         if not isinstance(other, SlidingWindowMultiplier):
             return False
-        return self.formulas == other.formulas and self.short_circuit == other.short_circuit and self.width == other.width and self.recoding_direction == other.recoding_direction and self.accumulation_order == other.accumulation_order
+        return (
+            self.formulas == other.formulas
+            and self.short_circuit == other.short_circuit
+            and self.width == other.width
+            and self.recoding_direction == other.recoding_direction
+            and self.accumulation_order == other.accumulation_order
+        )
 
     def __repr__(self):
         return f"{self.__class__.__name__}({', '.join(map(str, self.formulas.values()))}, short_circuit={self.short_circuit}, width={self.width}, recoding_direction={self.recoding_direction.name}, accumulation_order={self.accumulation_order.name})"
@@ -112,16 +129,20 @@ class FixedWindowLTRMultiplier(AccumulatorMultiplier, ScalarMultiplier):
     _points: MutableMapping[int, Point]
 
     def __init__(
-            self,
-            add: AdditionFormula,
-            dbl: DoublingFormula,
-            m: int,
-            scl: Optional[ScalingFormula] = None,
-            accumulation_order: AccumulationOrder = AccumulationOrder.PeqPR,
-            short_circuit: bool = True,
+        self,
+        add: AdditionFormula,
+        dbl: DoublingFormula,
+        m: int,
+        scl: Optional[ScalingFormula] = None,
+        accumulation_order: AccumulationOrder = AccumulationOrder.PeqPR,
+        short_circuit: bool = True,
     ):
         super().__init__(
-            short_circuit=short_circuit, accumulation_order=accumulation_order, add=add, dbl=dbl, scl=scl
+            short_circuit=short_circuit,
+            accumulation_order=accumulation_order,
+            add=add,
+            dbl=dbl,
+            scl=scl,
         )
         if m < 2:
             raise ValueError("Invalid base.")
@@ -134,7 +155,12 @@ class FixedWindowLTRMultiplier(AccumulatorMultiplier, ScalarMultiplier):
     def __eq__(self, other):
         if not isinstance(other, FixedWindowLTRMultiplier):
             return False
-        return self.formulas == other.formulas and self.short_circuit == other.short_circuit and self.m == other.m and self.accumulation_order == other.accumulation_order
+        return (
+            self.formulas == other.formulas
+            and self.short_circuit == other.short_circuit
+            and self.m == other.m
+            and self.accumulation_order == other.accumulation_order
+        )
 
     def __repr__(self):
         return f"{self.__class__.__name__}({', '.join(map(str, self.formulas.values()))}, short_circuit={self.short_circuit}, m={self.m}, accumulation_order={self.accumulation_order.name})"
@@ -180,3 +206,103 @@ class FixedWindowLTRMultiplier(AccumulatorMultiplier, ScalarMultiplier):
             if "scl" in self.formulas:
                 q = self._scl(q)
             return action.exit(q)
+
+
+@public
+class WindowBoothMultiplier(AccumulatorMultiplier, ScalarMultiplier):
+    """
+
+    :param short_circuit: Whether the use of formulas will be guarded by short-circuit on inputs
+                          of the point at infinity.
+    :param width: The width of the window.
+    :param accumulation_order: The order of accumulation of points.
+    :param precompute_negation: Whether to precompute the negation of the precomputed points as well.
+                                It is computed on the fly otherwise.
+    """
+
+    requires = {AdditionFormula, DoublingFormula, NegationFormula}
+    optionals = {ScalingFormula}
+    _points: MutableMapping[int, Point]
+    _points_neg: MutableMapping[int, Point]
+    precompute_negation: bool = False
+    """Whether to precompute the negation of the precomputed points as well."""
+    width: int
+    """The width of the window."""
+
+    def __init__(
+        self,
+        add: AdditionFormula,
+        dbl: DoublingFormula,
+        neg: NegationFormula,
+        width: int,
+        scl: Optional[ScalingFormula] = None,
+        accumulation_order: AccumulationOrder = AccumulationOrder.PeqPR,
+        precompute_negation: bool = False,
+        short_circuit: bool = True,
+    ):
+        super().__init__(
+            short_circuit=short_circuit,
+            accumulation_order=accumulation_order,
+            add=add,
+            dbl=dbl,
+            neg=neg,
+            scl=scl,
+        )
+        self.width = width
+        self.precompute_negation = precompute_negation
+
+    def __hash__(self):
+        return id(self)
+
+    def __eq__(self, other):
+        if not isinstance(other, WindowBoothMultiplier):
+            return False
+        return (
+            self.formulas == other.formulas
+            and self.short_circuit == other.short_circuit
+            and self.width == other.width
+            and self.precompute_negation == other.precompute_negation
+            and self.accumulation_order == other.accumulation_order
+        )
+
+    def __repr__(self):
+        return f"{self.__class__.__name__}({', '.join(map(str, self.formulas.values()))}, short_circuit={self.short_circuit}, width={self.width}, precompute_negation={self.precompute_negation}, accumulation_order={self.accumulation_order.name})"
+
+    def init(self, params: DomainParameters, point: Point):
+        with PrecomputationAction(params, point):
+            super().init(params, point)
+            double_point = self._dbl(point)
+            self._points = {1: point, 2: double_point}
+            if self.precompute_negation:
+                self._points_neg = {1: self._neg(point), 2: self._neg(double_point)}
+            current_point = double_point
+            for i in range(3, 2 ** (self.width - 1) + 1):
+                current_point = self._add(current_point, point)
+                self._points[i] = current_point
+                if self.precompute_negation:
+                    self._points_neg[i] = self._neg(current_point)
+
+    def multiply(self, scalar: int) -> Point:
+        if not self._initialized:
+            raise ValueError("ScalarMultiplier not initialized.")
+        with ScalarMultiplicationAction(self._point, scalar) as action:
+            if scalar == 0:
+                return action.exit(copy(self._params.curve.neutral))
+            scalar_booth = booth_window(
+                scalar, self.width, self._params.order.bit_length()
+            )
+            q = copy(self._params.curve.neutral)
+            for val in scalar_booth:
+                for _ in range(self.width):
+                    q = self._dbl(q)
+                if val > 0:
+                    q = self._accumulate(q, self._points[val])
+                elif val < 0:
+                    if self.precompute_negation:
+                        neg = self._points_neg[-val]
+                    else:
+                        neg = self._neg(self._points[-val])
+                    q = self._accumulate(q, neg)
+            if "scl" in self.formulas:
+                q = self._scl(q)
+            return action.exit(q)
diff --git a/pyecsca/ec/point.py b/pyecsca/ec/point.py
index 280d746..bdbba5e 100644
--- a/pyecsca/ec/point.py
+++ b/pyecsca/ec/point.py
@@ -140,7 +140,7 @@ class Point:
             if randomized:
                 lmbd = Mod.random(curve.prime)
                 for var, value in result.items():
-                    result[var] = value * lmbd**coordinate_model.homogweights[var]
+                    result[var] = value * (lmbd**coordinate_model.homogweights[var])
             return action.exit(Point(coordinate_model, **result))
 
     def equals_affine(self, other: "Point") -> bool:
diff --git a/pyecsca/ec/scalar.py b/pyecsca/ec/scalar.py
index 3fadc00..af5a6ab 100644
--- a/pyecsca/ec/scalar.py
+++ b/pyecsca/ec/scalar.py
@@ -1,5 +1,5 @@
 """Provides functions for computing various scalar representations (like NAF, or different bases)."""
-from typing import List
+from typing import List, Tuple, Literal
 from itertools import dropwhile
 from public import public
 
@@ -11,7 +11,7 @@ def convert_base(i: int, base: int) -> List[int]:
 
     :param i: The scalar.
     :param base: The base.
-    :return: The resulting digit list.
+    :return: The resulting digit list (little-endian).
     """
     if i == 0:
         return [0]
@@ -131,3 +131,62 @@ def naf(k: int) -> List[int]:
     :return: The NAF.
     """
     return wnaf(k, 2)
+
+
+@public
+def booth(k: int) -> List[int]:
+    """
+    Original Booth binary recoding, from [B51]_.
+
+    :param k: The scalar to recode.
+    :return: The recoded list of digits (0, 1, -1), little-endian.
+    """
+    res = []
+    for i in range(k.bit_length()):
+        a_i = (k >> i) & 1
+        b_i = (k >> (i + 1)) & 1
+        res.append(a_i - b_i)
+    res.insert(0, -(k & 1))
+    return res
+
+
+@public
+def booth_word(digit: int, w: int) -> int:
+    """
+    Modified Booth recoding, from [M61]_ and BoringSSL NIST impl.
+
+    Needs `w+1` bits of scalar in digit, but the one bit is overlapping (window size is `w`).
+
+    :param digit:
+    :param w:
+    :return:
+    """
+    if digit.bit_length() > (w + 1):
+        raise ValueError("Invalid digit, cannot be larger than w + 1 bits.")
+    s = ~((digit >> w) - 1)
+    d = (1 << (w + 1)) - digit - 1
+    d = (d & s) | (digit & ~s)
+    d = (d >> 1) + (d & 1)
+
+    return -d if s else d
+
+
+@public
+def booth_window(k: int, w: int, blen: int) -> List[int]:
+    """
+    Recode a whole scalar using Booth recoding as in BoringSSL.
+
+    :param k: The scalar.
+    :param w: The window size.
+    :param blen: The bit-length of the group.
+    :return: The big-endian recoding
+    """
+    mask = (1 << (w + 1)) - 1
+    res = []
+    for i in range(blen + (w - (blen % w) - 1), -1, -w):
+        if i >= w:
+            d = (k >> (i - w)) & mask
+        else:
+            d = (k << (w - i)) & mask
+        res.append(booth_word(d, w))
+    return res
diff --git a/pyecsca/sca/re/rpa.py b/pyecsca/sca/re/rpa.py
index db88dae..b1661e6 100644
--- a/pyecsca/sca/re/rpa.py
+++ b/pyecsca/sca/re/rpa.py
@@ -34,8 +34,11 @@ class MultipleContext(Context):
     """Context that traces the multiples of points computed."""
 
     base: Optional[Point]
+    """The base point that all the multiples are counted from."""
     points: MutableMapping[Point, int]
+    """The mapping of points to the multiples they represent (e.g., base -> 1)."""
     parents: MutableMapping[Point, List[Point]]
+    """The mapping of points to their parent they were computed from."""
     inside: bool
 
     def __init__(self):
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/pyecsca/sca/target/emulator.py b/pyecsca/sca/target/leakage.py
index bca27cf..8a8acd6 100644
--- a/pyecsca/sca/target/emulator.py
+++ b/pyecsca/sca/target/leakage.py
@@ -18,7 +18,7 @@ import numpy as np
 
 
 @public
-class EmulatorTarget(Target):
+class LeakageTarget(Target):
 
     model: CurveModel
     coords: CoordinateModel
@@ -48,7 +48,7 @@ class EmulatorTarget(Target):
         context.actions.walk(callback)
         return Trace(np.array(temp_trace))
 
-    def emulate_scalar_mult_traces(self, num_of_traces: int, scalar: int) -> Tuple[list[Point], list[Trace]]:
+    def simulate_scalar_mult_traces(self, num_of_traces: int, scalar: int) -> Tuple[list[Point], list[Trace]]:
         if self.params is None:
             raise ValueError("Missing DomainParameters")
         points = [self.params.curve.affine_random().to_model(self.coords, self.params.curve) for _ in range(num_of_traces)]
@@ -58,7 +58,7 @@ class EmulatorTarget(Target):
             traces.append(trace)
         return points, traces
 
-    def emulate_ecdh_traces(self, num_of_traces: int) -> Tuple[list[Point], list[Trace]]:
+    def simulate_ecdh_traces(self, num_of_traces: int) -> Tuple[list[Point], list[Trace]]:
         if self.params is None:
             raise ValueError("Missing DomainParameters")
         other_pubs = [self.params.curve.affine_random().to_model(self.coords, self.params.curve) for _ in range(num_of_traces)]
diff --git a/pyproject.toml b/pyproject.toml
index 52b03a1..48ad351 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -85,3 +85,6 @@ filterwarnings = [
 
 [tool.mypy]
 plugins = "numpy.typing.mypy_plugin"
+
+[tool.interrogate]
+exclude = ["pyecsca/ec/formula/fliparoo.py", "pyecsca/ec/formula/graph.py", "pyecsca/ec/formula/partitions.py", "pyecsca/ec/formula/switch_sign.py", "pyecsca/ec/std/.github/"]
diff --git a/test/ec/test_configuration.py b/test/ec/test_configuration.py
index 54cc53a..15c9471 100644
--- a/test/ec/test_configuration.py
+++ b/test/ec/test_configuration.py
@@ -33,7 +33,7 @@ def test_weierstrass_projective(base_independents):
     coords = model.coordinates["projective"]
     configs = list(all_configurations(model=model, coords=coords, **base_independents))
     assert len(set(map(lambda cfg: cfg.scalarmult, configs))) == len(configs)
-    assert len(configs) == 11360
+    assert len(configs) == 12640
 
 
 def test_mult_class(base_independents):
diff --git a/test/ec/test_formula.py b/test/ec/test_formula.py
index bb3d123..3f8d45c 100644
--- a/test/ec/test_formula.py
+++ b/test/ec/test_formula.py
@@ -1,13 +1,39 @@
 import pickle
+from operator import itemgetter
+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.expand import expand_formula_list
+from pyecsca.ec.formula.fliparoo import generate_fliparood_formulas
+from pyecsca.ec.formula.graph import rename_ivs
+from pyecsca.ec.formula.metrics import (
+    formula_similarity,
+    formula_similarity_abs,
+    formula_similarity_fuzz,
+)
+from pyecsca.ec.formula.partitions import (
+    reduce_all_adds,
+    expand_all_muls,
+    expand_all_nopower2_muls,
+    generate_partitioned_formulas,
+)
+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.efd import (
+    AdditionEFDFormula,
+    DoublingEFDFormula,
+    LadderEFDFormula,
+    EFDFormula,
+)
+from pyecsca.ec.formula import AdditionFormula, DoublingFormula, LadderFormula
 
 
 @pytest.fixture()
@@ -62,39 +88,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 +125,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 +138,371 @@ def test_symbolic(secp128r1, dbl):
 
 def test_pickle(add, dbl):
     assert add == pickle.loads(pickle.dumps(add))
+
+
+def test_formula_similarity(secp128r1):
+    add_bl = secp128r1.curve.coordinate_model.formulas["add-2007-bl"]
+    add_rcb = secp128r1.curve.coordinate_model.formulas["add-2015-rcb"]
+    out = formula_similarity(add_bl, add_rcb)
+    assert out["output"] == 0
+    assert out["ivs"] < 0.5
+    out_abs = formula_similarity_abs(add_bl, add_rcb)
+    assert out_abs["output"] == 0
+    assert out_abs["ivs"] < 0.5
+    out_fuzz = formula_similarity_fuzz(add_bl, add_rcb, secp128r1.curve, samples=100)
+    assert out_fuzz["output"] == 0
+    assert out_fuzz["ivs"] < 0.5
+    out_same = formula_similarity(add_bl, add_bl)
+    assert out_same["output"] == 1
+    assert out_same["ivs"] == 1
+
+
+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,
+    ],
+    [
+        "add-bearssl-v06",
+        ShortWeierstrassModel,
+        "jacobian",
+        ("secg", "secp256r1"),
+        AdditionEFDFormula,
+    ],
+    [
+        "dbl-bearssl-v06",
+        ShortWeierstrassModel,
+        "jacobian",
+        ("secg", "secp256r1"),
+        DoublingEFDFormula,
+    ],
+    [
+        "add-libgcrypt-v1102",
+        ShortWeierstrassModel,
+        "jacobian",
+        ("secg", "secp256r1"),
+        AdditionEFDFormula,
+    ],
+    [
+        "dbl-libgcrypt-v1102",
+        ShortWeierstrassModel,
+        "jacobian",
+        ("secg", "secp256r1"),
+        DoublingEFDFormula,
+    ],
+    [
+        "ladd-go-1214",
+        MontgomeryModel,
+        "xz",
+        ("other", "Curve25519"),
+        LadderEFDFormula,
+    ],
+    [
+        "add-gecc-322",
+        ShortWeierstrassModel,
+        "jacobian-3",
+        ("secg", "secp256r1"),
+        AdditionEFDFormula,
+    ],
+    [
+        "dbl-gecc-321",
+        ShortWeierstrassModel,
+        "jacobian-3",
+        ("secg", "secp256r1"),
+        DoublingEFDFormula,
+    ],
+    [
+        "ladd-boringssl-x25519",
+        MontgomeryModel,
+        "xz",
+        ("other", "Curve25519"),
+        LadderEFDFormula,
+    ],
+    [
+        "dbl-ipp-x25519",
+        MontgomeryModel,
+        "xz",
+        ("other", "Curve25519"),
+        DoublingEFDFormula,
+    ],
+    [
+        "ladd-botan-x25519",
+        MontgomeryModel,
+        "xz",
+        ("other", "Curve25519"),
+        LadderEFDFormula,
+    ],
+]
+
+
+@pytest.fixture(params=LIBRARY_FORMULAS, ids=list(map(itemgetter(0), LIBRARY_FORMULAS)))
+def library_formula_params(request) -> Tuple[EFDFormula, 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 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_partition_formula_single(library_formula_params):
+    formula, params = library_formula_params
+    try:
+        next(iter(generate_partitioned_formulas(formula)))
+    except StopIteration:
+        pass
+
+
+@pytest.mark.slow
+def test_partition_formula(library_formula_params):
+    formula, params = library_formula_params
+    for partitioned in generate_partitioned_formulas(formula):
+        do_test_formula(partitioned, 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):
+    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
+                )
+
+
+def test_formula_correctness(library_formula_params):
+    formula, params = library_formula_params
+    do_test_formula(formula, params)
+
+
+def test_formula_expand(add):
+    res = expand_formula_list([add])
+    assert len(res) > 1
diff --git a/test/ec/test_mult.py b/test/ec/test_mult.py
index 20e1c95..d5e3146 100644
--- a/test/ec/test_mult.py
+++ b/test/ec/test_mult.py
@@ -20,6 +20,7 @@ from pyecsca.ec.mult import (
     SlidingWindowMultiplier,
     BGMWMultiplier,
     CombMultiplier,
+    WindowBoothMultiplier,
 )
 from pyecsca.ec.mult.fixed import FullPrecompMultiplier
 from pyecsca.ec.point import InfinityPoint, Point
@@ -36,12 +37,9 @@ def assert_pt_equality(one: Point, other: Point, scale):
         assert one.equals(other)
 
 
-def do_basic_test(
-        mult_class, params, base, add, dbl, scale, neg=None, **kwargs
-):
+def do_basic_test(mult_class, params, base, add, dbl, scale, neg=None, **kwargs):
     mult = mult_class(
-        *get_formulas(params.curve.coordinate_model, add, dbl, neg, scale),
-        **kwargs
+        *get_formulas(params.curve.coordinate_model, add, dbl, neg, scale), **kwargs
     )
     mult.init(params, base)
     res = mult.multiply(314)
@@ -59,26 +57,28 @@ def do_basic_test(
     return res
 
 
-@pytest.mark.parametrize("add,dbl,scale",
-                         [
-                             ("add-1998-cmo-2", "dbl-1998-cmo-2", "z"),
-                             ("add-2015-rcb", "dbl-2015-rcb", None),
-                             ("add-1998-cmo-2", "dbl-1998-cmo-2", None),
-                         ])
+@pytest.mark.parametrize(
+    "add,dbl,scale",
+    [
+        ("add-1998-cmo-2", "dbl-1998-cmo-2", "z"),
+        ("add-2015-rcb", "dbl-2015-rcb", None),
+        ("add-1998-cmo-2", "dbl-1998-cmo-2", None),
+    ],
+)
 def test_rtl(secp128r1, add, dbl, scale):
     do_basic_test(RTLMultiplier, secp128r1, secp128r1.generator, add, dbl, scale)
 
 
-@pytest.mark.parametrize("add,dbl,scale",
-                         [
-                             ("add-1998-cmo-2", "dbl-1998-cmo-2", "z"),
-                             ("add-2015-rcb", "dbl-2015-rcb", None),
-                             ("add-1998-cmo-2", "dbl-1998-cmo-2", None),
-                         ])
+@pytest.mark.parametrize(
+    "add,dbl,scale",
+    [
+        ("add-1998-cmo-2", "dbl-1998-cmo-2", "z"),
+        ("add-2015-rcb", "dbl-2015-rcb", None),
+        ("add-1998-cmo-2", "dbl-1998-cmo-2", None),
+    ],
+)
 def test_ltr(secp128r1, add, dbl, scale):
-    a = do_basic_test(
-        LTRMultiplier, secp128r1, secp128r1.generator, add, dbl, scale
-    )
+    a = do_basic_test(LTRMultiplier, secp128r1, secp128r1.generator, add, dbl, scale)
     b = do_basic_test(
         LTRMultiplier, secp128r1, secp128r1.generator, add, dbl, scale, always=True
     )
@@ -100,22 +100,35 @@ def test_ltr(secp128r1, add, dbl, scale):
     assert_pt_equality(c, d, scale)
 
 
-@pytest.mark.parametrize("add,dbl,scale",
-                         [
-                             ("add-1998-cmo-2", "dbl-1998-cmo-2", "z"),
-                             ("add-2015-rcb", "dbl-2015-rcb", None),
-                             ("add-1998-cmo-2", "dbl-1998-cmo-2", None),
-                         ])
+@pytest.mark.parametrize(
+    "add,dbl,scale",
+    [
+        ("add-1998-cmo-2", "dbl-1998-cmo-2", "z"),
+        ("add-2015-rcb", "dbl-2015-rcb", None),
+        ("add-1998-cmo-2", "dbl-1998-cmo-2", None),
+    ],
+)
 def test_doubleandadd(secp128r1, add, dbl, scale):
     a = do_basic_test(
         DoubleAndAddMultiplier, secp128r1, secp128r1.generator, add, dbl, scale
     )
     b = do_basic_test(
-        DoubleAndAddMultiplier, secp128r1, secp128r1.generator, add, dbl, scale, direction=ProcessingDirection.RTL
+        DoubleAndAddMultiplier,
+        secp128r1,
+        secp128r1.generator,
+        add,
+        dbl,
+        scale,
+        direction=ProcessingDirection.RTL,
     )
     c = do_basic_test(
-        DoubleAndAddMultiplier, secp128r1, secp128r1.generator, add, dbl, scale,
-        accumulation_order=AccumulationOrder.PeqPR
+        DoubleAndAddMultiplier,
+        secp128r1,
+        secp128r1.generator,
+        add,
+        dbl,
+        scale,
+        accumulation_order=AccumulationOrder.PeqPR,
     )
     d = do_basic_test(
         DoubleAndAddMultiplier,
@@ -132,13 +145,14 @@ def test_doubleandadd(secp128r1, add, dbl, scale):
     assert_pt_equality(c, d, scale)
 
 
-@pytest.mark.parametrize("add,dbl,scale",
-                         [
-                             ("add-1998-cmo-2", "dbl-1998-cmo-2", "z"),
-                             ("add-2015-rcb", "dbl-2015-rcb", None),
-                             ("add-1998-cmo-2", "dbl-1998-cmo-2", None),
-                         ]
-                         )
+@pytest.mark.parametrize(
+    "add,dbl,scale",
+    [
+        ("add-1998-cmo-2", "dbl-1998-cmo-2", "z"),
+        ("add-2015-rcb", "dbl-2015-rcb", None),
+        ("add-1998-cmo-2", "dbl-1998-cmo-2", None),
+    ],
+)
 def test_coron(secp128r1, add, dbl, scale):
     do_basic_test(CoronMultiplier, secp128r1, secp128r1.generator, add, dbl, scale)
 
@@ -164,27 +178,31 @@ def test_ladder(curve25519):
     assert_pt_equality(a, b, True)
 
 
-@pytest.mark.parametrize("add,dbl,scale",
-                         [
-                             ("add-1998-cmo-2", "dbl-1998-cmo-2", "z"),
-                             ("add-2015-rcb", "dbl-2015-rcb", None),
-                             ("add-1998-cmo-2", "dbl-1998-cmo-2", None),
-                         ])
+@pytest.mark.parametrize(
+    "add,dbl,scale",
+    [
+        ("add-1998-cmo-2", "dbl-1998-cmo-2", "z"),
+        ("add-2015-rcb", "dbl-2015-rcb", None),
+        ("add-1998-cmo-2", "dbl-1998-cmo-2", None),
+    ],
+)
 def test_simple_ladder(secp128r1, add, dbl, scale):
     do_basic_test(
         SimpleLadderMultiplier, secp128r1, secp128r1.generator, add, dbl, scale
     )
 
 
-@pytest.mark.parametrize("num,complete",
-                         [
-                             (15, True),
-                             (15, False),
-                             (2355498743, True),
-                             (2355498743, False),
-                             (325385790209017329644351321912443757746, True),
-                             (325385790209017329644351321912443757746, False),
-                         ])
+@pytest.mark.parametrize(
+    "num,complete",
+    [
+        (15, True),
+        (15, False),
+        (2355498743, True),
+        (2355498743, False),
+        (325385790209017329644351321912443757746, True),
+        (325385790209017329644351321912443757746, False),
+    ],
+)
 def test_ladder_differential(curve25519, num, complete):
     ladder = LadderMultiplier(
         curve25519.curve.coordinate_model.formulas["ladd-1987-m"],
@@ -206,27 +224,31 @@ def test_ladder_differential(curve25519, num, complete):
     assert InfinityPoint(curve25519.curve.coordinate_model) == differential.multiply(0)
 
 
-@pytest.mark.parametrize("add,dbl,neg,scale",
-                         [
-                             ("add-1998-cmo-2", "dbl-1998-cmo-2", "neg", "z"),
-                             ("add-2015-rcb", "dbl-2015-rcb", "neg", None),
-                             ("add-1998-cmo-2", "dbl-1998-cmo-2", "neg", None),
-                         ])
+@pytest.mark.parametrize(
+    "add,dbl,neg,scale",
+    [
+        ("add-1998-cmo-2", "dbl-1998-cmo-2", "neg", "z"),
+        ("add-2015-rcb", "dbl-2015-rcb", "neg", None),
+        ("add-1998-cmo-2", "dbl-1998-cmo-2", "neg", None),
+    ],
+)
 def test_binary_naf(secp128r1, add, dbl, neg, scale):
     do_basic_test(
         BinaryNAFMultiplier, secp128r1, secp128r1.generator, add, dbl, scale, neg
     )
 
 
-@pytest.mark.parametrize("add,dbl,neg,width,scale",
-                         [
-                             ("add-1998-cmo-2", "dbl-1998-cmo-2", "neg", 3, "z"),
-                             ("add-1998-cmo-2", "dbl-1998-cmo-2", "neg", 3, None),
-                             ("add-2015-rcb", "dbl-2015-rcb", "neg", 3, None),
-                             ("add-1998-cmo-2", "dbl-1998-cmo-2", "neg", 5, "z"),
-                             ("add-1998-cmo-2", "dbl-1998-cmo-2", "neg", 5, None),
-                             ("add-2015-rcb", "dbl-2015-rcb", "neg", 5, None),
-                         ])
+@pytest.mark.parametrize(
+    "add,dbl,neg,width,scale",
+    [
+        ("add-1998-cmo-2", "dbl-1998-cmo-2", "neg", 3, "z"),
+        ("add-1998-cmo-2", "dbl-1998-cmo-2", "neg", 3, None),
+        ("add-2015-rcb", "dbl-2015-rcb", "neg", 3, None),
+        ("add-1998-cmo-2", "dbl-1998-cmo-2", "neg", 5, "z"),
+        ("add-1998-cmo-2", "dbl-1998-cmo-2", "neg", 5, None),
+        ("add-2015-rcb", "dbl-2015-rcb", "neg", 5, None),
+    ],
+)
 def test_window_naf(secp128r1, add, dbl, neg, width, scale):
     formulas = get_formulas(secp128r1.curve.coordinate_model, add, dbl, neg, scale)
     mult = WindowNAFMultiplier(*formulas[:3], width, *formulas[3:])
@@ -247,12 +269,14 @@ def test_window_naf(secp128r1, add, dbl, neg, width, scale):
     assert_pt_equality(res_precompute, res, scale)
 
 
-@pytest.mark.parametrize("add,dbl,width,scale",
-                         [
-                             ("add-1998-cmo-2", "dbl-1998-cmo-2", 5, "z"),
-                             ("add-2015-rcb", "dbl-2015-rcb", 5, None),
-                             ("add-1998-cmo-2", "dbl-1998-cmo-2", 5, None),
-                         ])
+@pytest.mark.parametrize(
+    "add,dbl,width,scale",
+    [
+        ("add-1998-cmo-2", "dbl-1998-cmo-2", 5, "z"),
+        ("add-2015-rcb", "dbl-2015-rcb", 5, None),
+        ("add-1998-cmo-2", "dbl-1998-cmo-2", 5, None),
+    ],
+)
 def test_fixed_window(secp128r1, add, dbl, width, scale):
     formulas = get_formulas(secp128r1.curve.coordinate_model, add, dbl, scale)
     mult = FixedWindowLTRMultiplier(*formulas[:2], width, *formulas[2:])
@@ -266,6 +290,27 @@ def test_fixed_window(secp128r1, add, dbl, width, scale):
     assert InfinityPoint(secp128r1.curve.coordinate_model) == mult.multiply(0)
 
 
+@pytest.mark.parametrize(
+    "add,dbl,neg,width,scale",
+    [
+        ("add-1998-cmo-2", "dbl-1998-cmo-2", "neg", 5, "z"),
+        ("add-2015-rcb", "dbl-2015-rcb", "neg", 5, None),
+        ("add-1998-cmo-2", "dbl-1998-cmo-2", "neg", 5, None),
+    ],
+)
+def test_booth_window(secp128r1, add, dbl, neg, width, scale):
+    formulas = get_formulas(secp128r1.curve.coordinate_model, add, dbl, neg, scale)
+    mult = WindowBoothMultiplier(*formulas[:3], width, *formulas[3:])
+    mult.init(secp128r1, secp128r1.generator)
+    res = mult.multiply(157 * 789)
+    other = mult.multiply(157)
+    mult.init(secp128r1, other)
+    other = mult.multiply(789)
+    assert_pt_equality(res, other, scale)
+    mult.init(secp128r1, secp128r1.generator)
+    assert InfinityPoint(secp128r1.curve.coordinate_model) == mult.multiply(0)
+
+
 @pytest.fixture(params=["add-1998-cmo-2", "add-2015-rcb"])
 def add(secp128r1, request):
     return secp128r1.curve.coordinate_model.formulas[request.param]
@@ -276,53 +321,137 @@ def dbl(secp128r1, request):
     return secp128r1.curve.coordinate_model.formulas[request.param]
 
 
-@pytest.mark.parametrize("num", [10, 2355498743, 325385790209017329644351321912443757746])
+@pytest.mark.parametrize(
+    "num", [10, 2355498743, 325385790209017329644351321912443757746]
+)
 def test_basic_multipliers(secp128r1, num, add, dbl):
     neg = secp128r1.curve.coordinate_model.formulas["neg"]
     scale = secp128r1.curve.coordinate_model.formulas["z"]
 
-    ltr_options = {"always": (True, False),
-                   "complete": (True, False),
-                   "accumulation_order": tuple(AccumulationOrder)}
-    ltrs = [LTRMultiplier(add, dbl, scale, **dict(zip(ltr_options.keys(), combination))) for combination in product(*ltr_options.values())]
+    ltr_options = {
+        "always": (True, False),
+        "complete": (True, False),
+        "accumulation_order": tuple(AccumulationOrder),
+    }
+    ltrs = [
+        LTRMultiplier(add, dbl, scale, **dict(zip(ltr_options.keys(), combination)))
+        for combination in product(*ltr_options.values())
+    ]
     rtl_options = ltr_options
-    rtls = [RTLMultiplier(add, dbl, scale, **dict(zip(rtl_options.keys(), combination))) for combination in product(*rtl_options.values())]
-    bnaf_options = {"direction": tuple(ProcessingDirection),
-                    "accumulation_order": tuple(AccumulationOrder)}
-    bnafs = [BinaryNAFMultiplier(add, dbl, neg, scale, **dict(zip(bnaf_options.keys(), combination))) for combination in product(*bnaf_options.values())]
-    wnaf_options = {"precompute_negation": (True, False),
-                    "width": (3, 5),
-                    "accumulation_order": tuple(AccumulationOrder)}
-    wnafs = [WindowNAFMultiplier(add, dbl, neg, scl=scale, **dict(zip(wnaf_options.keys(), combination))) for combination in product(*wnaf_options.values())]
+    rtls = [
+        RTLMultiplier(add, dbl, scale, **dict(zip(rtl_options.keys(), combination)))
+        for combination in product(*rtl_options.values())
+    ]
+    bnaf_options = {
+        "direction": tuple(ProcessingDirection),
+        "accumulation_order": tuple(AccumulationOrder),
+    }
+    bnafs = [
+        BinaryNAFMultiplier(
+            add, dbl, neg, scale, **dict(zip(bnaf_options.keys(), combination))
+        )
+        for combination in product(*bnaf_options.values())
+    ]
+    wnaf_options = {
+        "precompute_negation": (True, False),
+        "width": (3, 5),
+        "accumulation_order": tuple(AccumulationOrder),
+    }
+    wnafs = [
+        WindowNAFMultiplier(
+            add, dbl, neg, scl=scale, **dict(zip(wnaf_options.keys(), combination))
+        )
+        for combination in product(*wnaf_options.values())
+    ]
+    booth_options = {
+        "precompute_negation": (True, False),
+        "width": (3, 5),
+        "accumulation_order": tuple(AccumulationOrder),
+    }
+    booths = [
+        WindowBoothMultiplier(
+            add, dbl, neg, scl=scale, **dict(zip(booth_options.keys(), combination))
+        )
+        for combination in product(*booth_options.values())
+    ]
     ladder_options = {"complete": (True, False)}
-    ladders = [SimpleLadderMultiplier(add, dbl, scale, **dict(zip(ladder_options.keys(), combination))) for combination in product(*ladder_options.values())]
-    fixed_options = {"m": (5, 8),
-                     "accumulation_order": tuple(AccumulationOrder)}
-    fixeds = [FixedWindowLTRMultiplier(add, dbl, scl=scale, **dict(zip(fixed_options.keys(), combination))) for combination in product(*fixed_options.values())]
-    sliding_options = {"width": (3, 5),
-                       "recoding_direction": tuple(ProcessingDirection),
-                       "accumulation_order": tuple(AccumulationOrder)}
-    slides = [SlidingWindowMultiplier(add, dbl, scl=scale, **dict(zip(sliding_options.keys(), combination))) for combination in product(*sliding_options.values())]
-    precomp_options = {"always": (True, False),
-                       "complete": (True, False),
-                       "direction": tuple(ProcessingDirection),
-                       "accumulation_order": tuple(AccumulationOrder)}
-    precomps = [FullPrecompMultiplier(add, dbl, scl=scale, **dict(zip(precomp_options.keys(), combination))) for combination in product(*precomp_options.values())]
-    bgmw_options = {"width": (3, 5),
-                    "direction": tuple(ProcessingDirection),
-                    "accumulation_order": tuple(AccumulationOrder)}
-    bgmws = [BGMWMultiplier(add, dbl, scl=scale, **dict(zip(bgmw_options.keys(), combination))) for combination in product(*bgmw_options.values())]
-    comb_options = {"width": (2, 3, 5),
-                    "accumulation_order": tuple(AccumulationOrder)}
-    combs = [CombMultiplier(add, dbl, scl=scale, **dict(zip(comb_options.keys(), combination))) for combination in product(*comb_options.values())]
+    ladders = [
+        SimpleLadderMultiplier(
+            add, dbl, scale, **dict(zip(ladder_options.keys(), combination))
+        )
+        for combination in product(*ladder_options.values())
+    ]
+    fixed_options = {"m": (5, 8), "accumulation_order": tuple(AccumulationOrder)}
+    fixeds = [
+        FixedWindowLTRMultiplier(
+            add, dbl, scl=scale, **dict(zip(fixed_options.keys(), combination))
+        )
+        for combination in product(*fixed_options.values())
+    ]
+    sliding_options = {
+        "width": (3, 5),
+        "recoding_direction": tuple(ProcessingDirection),
+        "accumulation_order": tuple(AccumulationOrder),
+    }
+    slides = [
+        SlidingWindowMultiplier(
+            add, dbl, scl=scale, **dict(zip(sliding_options.keys(), combination))
+        )
+        for combination in product(*sliding_options.values())
+    ]
+    precomp_options = {
+        "always": (True, False),
+        "complete": (True, False),
+        "direction": tuple(ProcessingDirection),
+        "accumulation_order": tuple(AccumulationOrder),
+    }
+    precomps = [
+        FullPrecompMultiplier(
+            add, dbl, scl=scale, **dict(zip(precomp_options.keys(), combination))
+        )
+        for combination in product(*precomp_options.values())
+    ]
+    bgmw_options = {
+        "width": (3, 5),
+        "direction": tuple(ProcessingDirection),
+        "accumulation_order": tuple(AccumulationOrder),
+    }
+    bgmws = [
+        BGMWMultiplier(
+            add, dbl, scl=scale, **dict(zip(bgmw_options.keys(), combination))
+        )
+        for combination in product(*bgmw_options.values())
+    ]
+    comb_options = {"width": (2, 3, 5), "accumulation_order": tuple(AccumulationOrder)}
+    combs = [
+        CombMultiplier(
+            add, dbl, scl=scale, **dict(zip(comb_options.keys(), combination))
+        )
+        for combination in product(*comb_options.values())
+    ]
 
-    mults: Sequence[ScalarMultiplier] = ltrs + rtls + bnafs + wnafs + [CoronMultiplier(add, dbl, scale)] + ladders + fixeds + slides + precomps + bgmws + combs
+    mults: Sequence[ScalarMultiplier] = (
+        ltrs
+        + rtls
+        + bnafs
+        + wnafs
+        + booths
+        + [CoronMultiplier(add, dbl, scale)]
+        + ladders
+        + fixeds
+        + slides
+        + precomps
+        + bgmws
+        + combs
+    )
     results = []
     for mult in mults:
         mult.init(secp128r1, secp128r1.generator)
         res = mult.multiply(num)
         if results:
-            assert res == results[-1], f"Points not equal {res} != {results[-1]} for mult = {mult}"
+            assert (
+                res == results[-1]
+            ), f"Points not equal {res} != {results[-1]} for mult = {mult}"
         results.append(res)
 
 
diff --git a/test/ec/test_scalar.py b/test/ec/test_scalar.py
index 493690a..9c6cb60 100644
--- a/test/ec/test_scalar.py
+++ b/test/ec/test_scalar.py
@@ -1,4 +1,13 @@
-from pyecsca.ec.scalar import convert_base, sliding_window_ltr, sliding_window_rtl, wnaf, naf
+from pyecsca.ec.scalar import (
+    convert_base,
+    sliding_window_ltr,
+    sliding_window_rtl,
+    wnaf,
+    naf,
+    booth,
+    booth_word,
+    booth_window,
+)
 
 
 def test_convert():
@@ -14,3 +23,18 @@ def test_sliding_window():
 def test_nafs():
     i = 0b1100110101001101011011
     assert naf(i) == wnaf(i, 2)
+
+
+def test_booth():
+    assert booth(0b101) == [-1, 1, -1, 1]
+    for i in range(2**6):
+        bw = booth_word(i, 5)
+        # verified with BoringSSL recoding impl. (ec_GFp_nistp_recode_scalar_bits)
+        if i <= 31:
+            assert bw == (i + 1) // 2
+        else:
+            assert bw == -((64 - i) // 2)
+    s = 0x12345678123456781234567812345678123456781234567812345678
+    bw = booth_window(s, 5, 224)
+    # verified with BoringSSL ec_GFp_nistp224_point_mul
+    assert bw == [1, 4, 13, 3, -10, 15, 0, 9, 3, 9, -10, -12, -8, 2, 9, -6, 5, 13, -2, 1, -14, 7, -15, 11, 8, -16, 5, -14, -12, 11, -6, -4, 1, 4, 13, 3, -10, 15, 0, 9, 3, 9, -10, -12, -8]
diff --git a/test/sca/test_rpa.py b/test/sca/test_rpa.py
index f09a1a1..39c281d 100644
--- a/test/sca/test_rpa.py
+++ b/test/sca/test_rpa.py
@@ -12,12 +12,24 @@ from pyecsca.ec.mult import (
     RTLMultiplier,
     BinaryNAFMultiplier,
     WindowNAFMultiplier,
-    SimpleLadderMultiplier, AccumulationOrder, ProcessingDirection, SlidingWindowMultiplier, FixedWindowLTRMultiplier,
-    FullPrecompMultiplier, BGMWMultiplier, CombMultiplier,
+    SimpleLadderMultiplier,
+    AccumulationOrder,
+    ProcessingDirection,
+    SlidingWindowMultiplier,
+    FixedWindowLTRMultiplier,
+    FullPrecompMultiplier,
+    BGMWMultiplier,
+    CombMultiplier,
+    WindowBoothMultiplier,
 )
 from pyecsca.ec.params import DomainParameters
 from pyecsca.ec.point import Point
-from pyecsca.sca.re.rpa import MultipleContext, rpa_point_0y, rpa_point_x0, rpa_distinguish
+from pyecsca.sca.re.rpa import (
+    MultipleContext,
+    rpa_point_0y,
+    rpa_point_x0,
+    rpa_distinguish,
+)
 
 
 @pytest.fixture()
@@ -47,18 +59,18 @@ def neg(coords):
 
 @pytest.fixture()
 def rpa_params(model, coords):
-    p = 0x85d265945a4f5681
-    a = Mod(0x7fc57b4110698bc0, p)
-    b = Mod(0x37113ea591b04527, p)
-    gx = Mod(0x80d2d78fddb97597, p)
-    gy = Mod(0x5586d818b7910930, p)
+    p = 0x85D265945A4F5681
+    a = Mod(0x7FC57B4110698BC0, p)
+    b = Mod(0x37113EA591B04527, p)
+    gx = Mod(0x80D2D78FDDB97597, p)
+    gy = Mod(0x5586D818B7910930, p)
     # (0x4880bcf620852a54, 0) RPA point
     # (0, 0x6bed3155c9ada064) RPA point
 
     infty = Point(coords, X=Mod(0, p), Y=Mod(1, p), Z=Mod(0, p))
     g = Point(coords, X=gx, Y=gy, Z=Mod(1, p))
     curve = EllipticCurve(model, coords, p, infty, dict(a=a, b=b))
-    return DomainParameters(curve, g, 0x85d265932d90785c, 1)
+    return DomainParameters(curve, g, 0x85D265932D90785C, 1)
 
 
 def test_x0_point(rpa_params):
@@ -80,25 +92,63 @@ def test_distinguish(secp128r1, add, dbl, neg):
         RTLMultiplier(add, dbl, None, False, AccumulationOrder.PeqPR, True),
         RTLMultiplier(add, dbl, None, True, AccumulationOrder.PeqPR, False),
         SimpleLadderMultiplier(add, dbl, None, True, True),
-        BinaryNAFMultiplier(add, dbl, neg, None, ProcessingDirection.LTR, AccumulationOrder.PeqPR, True),
-        WindowNAFMultiplier(add, dbl, neg, 3, None, AccumulationOrder.PeqPR, True, True),
-        WindowNAFMultiplier(add, dbl, neg, 4, None, AccumulationOrder.PeqPR, True, True),
+        BinaryNAFMultiplier(
+            add, dbl, neg, None, ProcessingDirection.LTR, AccumulationOrder.PeqPR, True
+        ),
+        WindowNAFMultiplier(
+            add, dbl, neg, 3, None, AccumulationOrder.PeqPR, True, True
+        ),
+        WindowNAFMultiplier(
+            add, dbl, neg, 4, None, AccumulationOrder.PeqPR, True, True
+        ),
+        WindowBoothMultiplier(
+            add, dbl, neg, 4, None, AccumulationOrder.PeqPR, True, True
+        ),
         # WindowNAFMultiplier(add, dbl, neg, 4, None, AccumulationOrder.PeqPR, False, True),
-        SlidingWindowMultiplier(add, dbl, 3, None, ProcessingDirection.LTR, AccumulationOrder.PeqPR, True),
-        SlidingWindowMultiplier(add, dbl, 5, None, ProcessingDirection.LTR, AccumulationOrder.PeqPR, True),
+        SlidingWindowMultiplier(
+            add, dbl, 3, None, ProcessingDirection.LTR, AccumulationOrder.PeqPR, True
+        ),
+        SlidingWindowMultiplier(
+            add, dbl, 5, None, ProcessingDirection.LTR, AccumulationOrder.PeqPR, True
+        ),
         FixedWindowLTRMultiplier(add, dbl, 4, None, AccumulationOrder.PeqPR, True),
         FixedWindowLTRMultiplier(add, dbl, 5, None, AccumulationOrder.PeqPR, True),
-        FullPrecompMultiplier(add, dbl, None, True, ProcessingDirection.LTR, AccumulationOrder.PeqPR, True, True),
-        FullPrecompMultiplier(add, dbl, None, False, ProcessingDirection.LTR, AccumulationOrder.PeqPR, True, True),
+        FullPrecompMultiplier(
+            add,
+            dbl,
+            None,
+            True,
+            ProcessingDirection.LTR,
+            AccumulationOrder.PeqPR,
+            True,
+            True,
+        ),
+        FullPrecompMultiplier(
+            add,
+            dbl,
+            None,
+            False,
+            ProcessingDirection.LTR,
+            AccumulationOrder.PeqPR,
+            True,
+            True,
+        ),
         # FullPrecompMultiplier(add, dbl, None, False, ProcessingDirection.RTL, AccumulationOrder.PeqPR, True, True),
-        BGMWMultiplier(add, dbl, 3, None, ProcessingDirection.LTR, AccumulationOrder.PeqPR, True),
-        BGMWMultiplier(add, dbl, 5, None, ProcessingDirection.LTR, AccumulationOrder.PeqPR, True),
+        BGMWMultiplier(
+            add, dbl, 3, None, ProcessingDirection.LTR, AccumulationOrder.PeqPR, True
+        ),
+        BGMWMultiplier(
+            add, dbl, 5, None, ProcessingDirection.LTR, AccumulationOrder.PeqPR, True
+        ),
         CombMultiplier(add, dbl, 3, None, AccumulationOrder.PeqPR, True),
-        CombMultiplier(add, dbl, 5, None, AccumulationOrder.PeqPR, True)
+        CombMultiplier(add, dbl, 5, None, AccumulationOrder.PeqPR, True),
     ]
     for real_mult in multipliers:
+
         def simulated_oracle(scalar, affine_point):
-            point = affine_point.to_model(secp128r1.curve.coordinate_model, secp128r1.curve)
+            point = affine_point.to_model(
+                secp128r1.curve.coordinate_model, secp128r1.curve
+            )
             with local(MultipleContext()) as ctx:
                 real_mult.init(secp128r1, point)
                 real_mult.multiply(scalar)
diff --git a/test/sca/test_structural.py b/test/sca/test_structural.py
deleted file mode 100644
index 970e4fc..0000000
--- a/test/sca/test_structural.py
+++ /dev/null
@@ -1,315 +0,0 @@
-import pytest
-from importlib_resources import files, as_file
-import test.data.formulas
-from pyecsca.ec.formula import (
-    AdditionEFDFormula,
-    LadderEFDFormula,
-    DoublingEFDFormula,
-    AdditionFormula,
-    DoublingFormula,
-    LadderFormula,
-)
-from pyecsca.ec.model import ShortWeierstrassModel, MontgomeryModel, TwistedEdwardsModel
-from pyecsca.ec.params import get_params
-from pyecsca.sca.re.structural import formula_similarity
-
-
-def test_formula_similarity(secp128r1):
-    add_bl = secp128r1.curve.coordinate_model.formulas["add-2007-bl"]
-    add_rcb = secp128r1.curve.coordinate_model.formulas["add-2015-rcb"]
-    out = formula_similarity(add_bl, add_rcb)
-    assert out["output"] == 0
-    assert out["ivs"] < 0.5
-    out_same = formula_similarity(add_bl, add_bl)
-    assert out_same["output"] == 1
-    assert out_same["ivs"] == 1
-
-
-@pytest.mark.parametrize(
-    "name,model,coords,param_spec,formula_type",
-    [
-        [
-            "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,
-        ],
-        [
-            "add-bearssl-v06",
-            ShortWeierstrassModel,
-            "jacobian",
-            ("secg", "secp256r1"),
-            AdditionEFDFormula,
-        ],
-        [
-            "dbl-bearssl-v06",
-            ShortWeierstrassModel,
-            "jacobian",
-            ("secg", "secp256r1"),
-            DoublingEFDFormula,
-        ],
-        [
-            "add-libgcrypt-v1102",
-            ShortWeierstrassModel,
-            "jacobian",
-            ("secg", "secp256r1"),
-            AdditionEFDFormula,
-        ],
-        [
-            "dbl-libgcrypt-v1102",
-            ShortWeierstrassModel,
-            "jacobian",
-            ("secg", "secp256r1"),
-            DoublingEFDFormula,
-        ],
-        [
-            "ladd-go-1214",
-            MontgomeryModel,
-            "xz",
-            ("other", "Curve25519"),
-            LadderEFDFormula,
-        ],
-        [
-            "add-gecc-322",
-            ShortWeierstrassModel,
-            "jacobian-3",
-            ("secg", "secp256r1"),
-            AdditionEFDFormula,
-        ],
-        [
-            "dbl-gecc-321",
-            ShortWeierstrassModel,
-            "jacobian-3",
-            ("secg", "secp256r1"),
-            DoublingEFDFormula,
-        ],
-        [
-            "ladd-boringssl-x25519",
-            MontgomeryModel,
-            "xz",
-            ("other", "Curve25519"),
-            LadderEFDFormula,
-        ],
-        [
-            "dbl-ipp-x25519",
-            MontgomeryModel,
-            "xz",
-            ("other", "Curve25519"),
-            DoublingEFDFormula,
-        ],
-        [
-            "ladd-botan-x25519",
-            MontgomeryModel,
-            "xz",
-            ("other", "Curve25519"),
-            LadderEFDFormula,
-        ],
-    ],
-)
-def test_formula_correctness(name, model, coords, param_spec, formula_type):
-    model = model()
-    coordinate_model = model.coordinates[coords]
-    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)
-    scale = coordinate_model.formulas.get("z", None)
-    if scale is None:
-        scale = coordinate_model.formulas.get("scale", None)
-    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
-                )