Implement multipliers.

26 files changed, 516 insertions, 71 deletions

diff --git a/pyecsca/codegen/bn/bn.c b/pyecsca/codegen/bn/bn.c
index dd22ca1..635eaae 100644
--- a/pyecsca/codegen/bn/bn.c
+++ b/pyecsca/codegen/bn/bn.c
@@ -1,4 +1,6 @@
 #include "bn.h"
+#include <string.h>
+#include <stdlib.h>
 
 bn_err bn_init(bn_t *bn) {
 	return mp_init(bn);
@@ -28,6 +30,7 @@ int bn_from_int(uint64_t value, bn_t *out) {
 void bn_to_binpad(const bn_t *one, uint8_t *data, size_t size) {
 	size_t ubin_size = mp_ubin_size(one);
 	size_t offset = size - ubin_size;
+	memset(data, 0, offset);
 	mp_to_ubin(one, data + offset, ubin_size, NULL);
 }
 
@@ -109,6 +112,18 @@ bool bn_eq(const bn_t *one, const bn_t *other) {
 	return mp_cmp_mag(one, other) == MP_EQ;
 }
 
+bool bn_is_0(const bn_t *one) {
+	return mp_cmp_d(one, 0) == MP_EQ;
+}
+
+bool bn_is_1(const bn_t *one) {
+	return mp_cmp_d(one, 1) == MP_EQ;
+}
+
+bn_sign bn_get_sign(const bn_t *one) {
+	return one->sign;
+}
+
 int bn_get_bit(const bn_t *bn, int which) {
 	int which_digit = which / MP_DIGIT_BIT;
 	int which_bit = which % MP_DIGIT_BIT;
@@ -120,4 +135,53 @@ int bn_get_bit(const bn_t *bn, int which) {
 
 int bn_bit_length(const bn_t *bn) {
 	return mp_count_bits(bn);
+}
+
+wnaf_t *bn_wnaf(const bn_t *bn, int w) {
+	if (w > 8 || w < 2) {
+		return NULL;
+	}
+	wnaf_t *result = malloc(sizeof(wnaf_t));
+	result->w = w;
+	result->length = bn_bit_length(bn) + 1;
+	result->data = calloc(result->length, sizeof(int8_t));
+
+	bn_t half_width;
+	bn_init(&half_width);
+	bn_from_int(1, &half_width);
+	bn_lsh(&half_width, w - 1, &half_width);
+	bn_t full_width;
+	bn_init(&full_width);
+	bn_from_int(1, &full_width);
+	bn_lsh(&full_width, w, &full_width);
+
+	bn_t k; bn_init(&k);
+	bn_copy(bn, &k);
+
+	bn_t val_mod; bn_init(&val_mod);
+
+	size_t i = 0;
+	while (!bn_is_0(&k) && !(bn_get_sign(&k) == MP_NEG)) {
+		if (bn_get_bit(&k, 0) == 1) {
+			bn_mod(&k, &full_width, &val_mod);
+			if (mp_cmp(&val_mod, &half_width) == MP_GT) {
+				mp_sub(&val_mod, &full_width, &val_mod);
+			}
+			int8_t val = (int8_t) mp_get_i32(&val_mod);
+			result->data[i++] = val;
+			mp_sub(&k, &val_mod, &k);
+		} else {
+			result->data[i++] = 0;
+		}
+		bn_rsh(&k, 1, &k);
+	}
+	bn_clear(&val_mod);
+	bn_clear(&half_width);
+	bn_clear(&full_width);
+	bn_clear(&k);
+	return result;
+}
+
+wnaf_t *bn_bnaf(const bn_t *bn) {
+	return bn_wnaf(bn, 2);
 }
\ No newline at end of file
diff --git a/pyecsca/codegen/bn/bn.h b/pyecsca/codegen/bn/bn.h
index 68ede2f..a5b1ba9 100644
--- a/pyecsca/codegen/bn/bn.h
+++ b/pyecsca/codegen/bn/bn.h
@@ -5,12 +5,19 @@
 
 #define bn_t mp_int
 #define bn_err mp_err
+#define bn_sign mp_sign
 
 typedef struct {
 	char name;
 	bn_t value;
 } named_bn_t;
 
+typedef struct {
+	int8_t *data;
+	size_t length;
+	int w;
+} wnaf_t;
+
 bn_err bn_init(bn_t *bn);
 void bn_copy(const bn_t *from, bn_t *to);
 void bn_clear(bn_t *bn);
@@ -45,8 +52,13 @@ void bn_lsh(const bn_t *one, int amount, bn_t *out);
 void bn_rsh(const bn_t *one, int amount, bn_t *out);
 
 bool bn_eq(const bn_t *one, const bn_t *other);
+bool bn_is_0(const bn_t *one);
+bool bn_is_1(const bn_t *one);
+bn_sign bn_get_sign(const bn_t *one);
 
 int bn_get_bit(const bn_t *bn, int which);
 int bn_bit_length(const bn_t *bn);
+wnaf_t *bn_wnaf(const bn_t *bn, int w);
+wnaf_t *bn_bnaf(const bn_t *bn);
 
 #endif //BN_H_
\ No newline at end of file
diff --git a/pyecsca/codegen/builder.py b/pyecsca/codegen/builder.py
index d233bdc..366ee18 100644
--- a/pyecsca/codegen/builder.py
+++ b/pyecsca/codegen/builder.py
@@ -1,9 +1,11 @@
 #!/usr/bin/env python3
 import os
+import re
 import shutil
 import subprocess
 import tempfile
 from ast import operator, Add, Sub, Mult, Div, Pow
+from copy import copy
 from dataclasses import dataclass
 from enum import Enum
 from os import path
@@ -18,6 +20,9 @@ from pyecsca.ec.formula import (Formula, AdditionFormula, DoublingFormula, Tripl
                                 LadderFormula)
 from pyecsca.ec.model import (CurveModel, ShortWeierstrassModel, MontgomeryModel, EdwardsModel,
                               TwistedEdwardsModel)
+from pyecsca.ec.mult import (ScalarMultiplier, LTRMultiplier, RTLMultiplier, CoronMultiplier,
+                             LadderMultiplier, SimpleLadderMultiplier, DifferentialLadderMultiplier,
+                             WindowNAFMultiplier, BinaryNAFMultiplier)
 from pyecsca.ec.op import CodeOp
 
 env = Environment(
@@ -38,11 +43,14 @@ def render_op(op: operator, result: str, left: str, right: str, mod: str) -> Opt
         return "bn_mod_sqr(&{}, &{}, &{});".format(left, mod, result)
     elif isinstance(op, Pow):
         return "bn_mod_pow(&{}, &{}, &{}, &{});".format(left, right, mod, result)
+    elif op is None:
+        return "bn_copy(&{}, &{});".format(left, result)
     else:
         print(op, result, left, right, mod)
 
 
 env.globals["render_op"] = render_op
+env.globals["isinstance"] = isinstance
 
 
 class EnumDefine(Enum):
@@ -58,6 +66,7 @@ class EnumDefine(Enum):
 
 
 class Platform(EnumDefine):
+    """Platform to build for."""
     HOST = "HOST"
     XMEGA = "CW308_XMEGA"
     STM32F0 = "CW308_STM32F0"
@@ -65,6 +74,7 @@ class Platform(EnumDefine):
 
 
 class Multiplication(EnumDefine):
+    """Base multiplication algorithm to use."""
     TOOM_COOK = "MUL_TOOM_COOK"
     KARATSUBA = "MUL_KARATSUBA"
     COMBA = "MUL_COMBA"
@@ -72,6 +82,7 @@ class Multiplication(EnumDefine):
 
 
 class Squaring(EnumDefine):
+    """Base squaring algorithm to use."""
     TOOM_COOK = "SQR_TOOM_COOK"
     KARATSUBA = "SQR_KARATSUBA"
     COMBA = "SQR_COMBA"
@@ -79,12 +90,14 @@ class Squaring(EnumDefine):
 
 
 class Reduction(EnumDefine):
+    """Modular reduction method used."""
     BARRETT = "RED_BARRETT"
     MONTGOMERY = "RED_MONTGOMERY"
     BASE = "RED_BASE"
 
 
 class HashType(EnumDefine):
+    """Hash algorithm used in ECDH and ECDSA."""
     NONE = "HASH_NONE"
     SHA1 = "HASH_SHA1"
     SHA224 = "HASH_SHA224"
@@ -93,12 +106,8 @@ class HashType(EnumDefine):
     SHA512 = "HASH_SHA512"
 
 
-class ScalarMultAlgo(EnumDefine):
-    NONE = "MULT_NONE"
-    DOUBLE_AND_ADD = "MULT_DOUBLE_AND_ADD"
-
-
 class RandomMod(EnumDefine):
+    """Method of sampling a uniform integer modulo order."""
     SAMPLE = "MOD_RAND_SAMPLE"
     REDUCE = "MOD_RAND_REDUCE"
 
@@ -107,7 +116,6 @@ class RandomMod(EnumDefine):
 class Configuration(object):
     platform: Platform
     hash_type: HashType
-    mult_algo: ScalarMultAlgo
     mod_rand: RandomMod
     mult: Multiplication  # TODO: Use this
     sqr: Squaring  # TODO: Use this
@@ -115,6 +123,7 @@ class Configuration(object):
     model: CurveModel
     coords: CoordinateModel
     formulas: List[Formula]
+    scalarmult: ScalarMultiplier
 
 
 def render_defs(model: CurveModel, coords: CoordinateModel) -> str:
@@ -203,7 +212,7 @@ def render_coords_impl(coords: CoordinateModel) -> str:
                                               to_affine_rets=returns, to_affine_frees=frees)
 
 
-def render_formula_impl(formula: Formula) -> str:
+def render_formula_impl(formula: Formula, short_circuit: bool = False) -> str:
     if isinstance(formula, AdditionFormula):
         tname = "formula_add.c"
     elif isinstance(formula, DoublingFormula):
@@ -236,18 +245,26 @@ def render_formula_impl(formula: Formula) -> str:
         renames[output] = "{}->{}".format(outputs[num], var)
     namespace = transform_ops(formula.code, formula.coordinate_model.curve_model.parameter_names,
                               formula.outputs, renames)
+    namespace["short_circuit"] = short_circuit
     return template.render(namespace)
 
 
+def render_scalarmult_impl(scalarmult: ScalarMultiplier) -> str:
+    return env.get_template("mult.c").render(scalarmult=scalarmult, LTRMultiplier=LTRMultiplier,
+                                             RTLMultiplier=RTLMultiplier, CoronMultiplier=CoronMultiplier,
+                                             LadderMultiplier=LadderMultiplier, SimpleLadderMultiplier=SimpleLadderMultiplier,
+                                             DifferentialLadderMultiplier=DifferentialLadderMultiplier,
+                                             BinaryNAFMultiplier=BinaryNAFMultiplier)
+
+
 def render_main(model: CurveModel, coords: CoordinateModel) -> str:
     return env.get_template("main.c").render(curve_variables=coords.variables,
                                              curve_parameters=model.parameter_names)
 
 
-def render_makefile(platform: Platform, hash_type: HashType, mult_algo: ScalarMultAlgo,
-                    mod_rand: RandomMod) -> str:
+def render_makefile(platform: Platform, hash_type: HashType, mod_rand: RandomMod) -> str:
     return env.get_template("Makefile").render(platform=str(platform), hash_type=str(hash_type),
-                                               mult_algo=str(mult_algo), mod_rand=str(mod_rand))
+                                               mod_rand=str(mod_rand))
 
 
 def save_render(dir: str, fname: str, rendered: str):
@@ -258,22 +275,22 @@ def save_render(dir: str, fname: str, rendered: str):
 def render(config: Configuration) -> Tuple[str, str]:
     temp = tempfile.mkdtemp()
     symlinks = ["asn1", "bn", "hal", "hash", "mult", "prng", "simpleserial", "tommath", "fat.h",
-                "point.h", "curve.h", "Makefile.inc"]
+                "point.h", "curve.h", "mult.h", "Makefile.inc"]
     for sym in symlinks:
         os.symlink(resource_filename("pyecsca.codegen", sym), path.join(temp, sym))
     gen_dir = path.join(temp, "gen")
     os.mkdir(gen_dir)
     save_render(temp, "Makefile",
-                render_makefile(config.platform, config.hash_type, config.mult_algo,
-                                config.mod_rand))
+                render_makefile(config.platform, config.hash_type, config.mod_rand))
     save_render(temp, "main.c", render_main(config.model, config.coords))
     save_render(gen_dir, "defs.h", render_defs(config.model, config.coords))
     point_render = render_coords_impl(config.coords)
     for formula in config.formulas:
         point_render += "\n"
-        point_render += render_formula_impl(formula)
+        point_render += render_formula_impl(formula, config.scalarmult.short_circuit)
     save_render(gen_dir, "point.c", point_render)
     save_render(gen_dir, "curve.c", render_curve_impl(config.model))
+    save_render(gen_dir, "mult.c", render_scalarmult_impl(config.scalarmult))
     return temp, "pyecsca-codegen-{}.elf".format(str(config.platform))
 
 
@@ -314,9 +331,60 @@ def get_formula(ctx: click.Context, param, value: Optional[Tuple[str]]) -> List[
             raise click.BadParameter(
                     "Formula '{}' is not a formula in '{}'.".format(formula, coords))
         result.append(coords.formulas[formula])
+    if len(set(formula.__class__ for formula in result)) != len(result):
+        raise click.BadParameter("Duplicate formula types.")
+    ctx.meta["formulas"] = copy(result)
     return result
 
 
+def get_multiplier(ctx: click.Context, param, value: Optional[str]) -> Optional[ScalarMultiplier]:
+    if value is None:
+        return None
+    res = re.match(
+        "(?P<name>[a-zA-Z\-]+)\((?P<args>([a-zA-Z_]+ *= *[a-zA-Z0-9]+, )*?([a-zA-Z_]+ *= *[a-zA-Z0-9]+)*)\)",
+        value)
+    if not res:
+        raise click.BadParameter("Couldn't parse multiplier spec: {}.".format(value))
+    name = res.group("name")
+    args = res.group("args")
+    if name in ("ltr", "LTRMultiplier"):
+        mult_class = LTRMultiplier
+    elif name in ("rtl", "RTLMultiplier"):
+        mult_class = RTLMultiplier
+    elif name in ("coron", "CoronMultiplier"):
+        mult_class = CoronMultiplier
+    elif name in ("ldr", "LadderMultiplier"):
+        mult_class = LadderMultiplier
+    elif name in ("simple-ldr", "SimpleLadderMultiplier"):
+        mult_class = SimpleLadderMultiplier
+    elif name in ("diff-ldr", "DifferentialLadderMultiplier"):
+        mult_class = DifferentialLadderMultiplier
+    elif name in ("naf", "bnaf", "BinaryNAFMultiplier"):
+        mult_class = BinaryNAFMultiplier
+    elif name in ("wnaf", "WindowNAFMultiplier"):
+        mult_class = WindowNAFMultiplier
+    else:
+        raise click.BadParameter("Unknown multiplier: {}.".format(name))
+    formulas = ctx.meta["formulas"]
+    classes = set(formula.__class__ for formula in formulas)
+    if not all(
+            any(issubclass(cls, required) for cls in classes) for required in mult_class.requires):
+        raise click.BadParameter(
+            "Multiplier {} requires formulas: {}, got {}.".format(mult_class.__name__,
+                                                                  mult_class.requires, classes))
+    kwargs = eval("dict(" + args + ")")
+    required = set(
+        filter(lambda formula: any(isinstance(formula, cls) for cls in mult_class.requires),
+               formulas))
+    optional = set(
+        filter(lambda formula: any(isinstance(formula, cls) for cls in mult_class.optionals),
+               formulas))
+    for formula in required.union(optional):
+        kwargs[formula.shortname] = formula
+    mult = mult_class(**kwargs)
+    return mult
+
+
 def wrap_enum(enum_class: Type[EnumDefine]):
     def callback(ctx, param, value):
         try:
@@ -344,10 +412,6 @@ def main():
               type=click.Choice(HashType.names()),
               callback=wrap_enum(HashType),
               help="The hash algorithm to use (in ECDH and ECDSA).")
-@click.option("--mult", envvar="MULT_ALGO", required=True,
-              type=click.Choice(ScalarMultAlgo.names()),
-              callback=wrap_enum(ScalarMultAlgo),
-              help="The scalar multiplication algorithm to use.")
 @click.option("--rand", envvar="MOD_RAND", default="SAMPLE", show_default=True,
               type=click.Choice(RandomMod.names()),
               callback=wrap_enum(RandomMod),
@@ -373,17 +437,22 @@ def main():
                 callback=get_coords)
 @click.argument("formulas", required=True, nargs=-1,
                 callback=get_formula)
+@click.argument("scalarmult", required=True,
+                callback=get_multiplier)
 @click.argument("outfile")
-def build(platform, hash, mult, rand, mul, sqr, red, strip, remove, model, coords, formulas, outfile):
+def build(platform, hash, rand, mul, sqr, red, strip, remove, model, coords, formulas, scalarmult,
+          outfile):
     """This command builds an ECC implementation.
 
     \b
     MODEL: The curve model to use.
     COORDS: The coordinate model to use.
     FORMULAS: The formulas to use.
+    MULT: The scalar multiplication algorithm to use.
     OUTFILE: The output binary file with the built impl.
     """
-    config = Configuration(platform, hash, mult, rand, mul, sqr, red, model, coords, formulas)
+
+    config = Configuration(platform, hash, rand, mul, sqr, red, model, coords, formulas, scalarmult)
     dir, file = render(config)
 
     full_path = path.join(dir, file)
diff --git a/pyecsca/codegen/client.py b/pyecsca/codegen/client.py
index 9d839fb..964487b 100644
--- a/pyecsca/codegen/client.py
+++ b/pyecsca/codegen/client.py
@@ -8,6 +8,7 @@ from pyecsca.ec.curve import EllipticCurve
 from pyecsca.ec.group import AbelianGroup
 from pyecsca.ec.mod import Mod
 from pyecsca.ec.model import ShortWeierstrassModel
+from pyecsca.ec.mult import LTRMultiplier
 from pyecsca.ec.point import Point
 
 
@@ -67,24 +68,71 @@ def cmd_set_curve(group: AbelianGroup) -> str:
     return "c" + hexlify(encode_data(None, data)).decode()
 
 
+def cmd_generate() -> str:
+    return "g"
+
+
+def cmd_set_privkey(privkey: int) -> str:
+    return "s" + hexlify(encode_data(None, {"s": encode_scalar(privkey)})).decode()
+
+
+def cmd_set_pubkey(pubkey: Point) -> str:
+    return "w" + hexlify(encode_data(None, {"w": encode_point(pubkey)})).decode()
+
+
+def cmd_scalar_mult(scalar: int) -> str:
+    return "m" + hexlify(encode_data(None, {"s": encode_scalar(scalar)})).decode()
+
+
+def cmd_ecdh(pubkey: Point) -> str:
+    return "e" + hexlify(encode_data(None, {"w": encode_point(pubkey)})).decode()
+
+
+def cmd_ecdsa_sign(data: bytes) -> str:
+    return "a" + hexlify(encode_data(None, {"d": data})).decode()
+
+
+def cmd_ecdsa_verify(data: bytes, sig: bytes) -> str:
+    return "v" + hexlify(encode_data(None, {"d": data, "s": sig})).decode()
+
+
 @click.command()
 @click.version_option()
 def main():
-    model = ShortWeierstrassModel()
-    coords = model.coordinates["projective"]
-    p = 0xffffffff00000001000000000000000000000000ffffffffffffffffffffffff
-    curve = EllipticCurve(model, coords,
-                          p,
-                          {"a": 0xffffffff00000001000000000000000000000000fffffffffffffffffffffffc,
-                           "b": 0x5ac635d8aa3a93e7b3ebbd55769886bc651d06b0cc53b0f63bce3c3e27d2604b})
-    affine_g = Point(AffineCoordinateModel(model),
-                     x=Mod(0x6b17d1f2e12c4247f8bce6e563a440f277037d812deb33a0f4a13945d898c296, p),
-                     y=Mod(0x4fe342e2fe1a7f9b8ee7eb4a7c0f9e162bce33576b315ececbb6406837bf51f5, p))
-    g = Point.from_affine(coords, affine_g)
-    neutral = Point(coords, X=Mod(0, p), Y=Mod(1, p), Z=Mod(0, p))
-    group = AbelianGroup(curve, g, neutral,
-                         0xffffffff00000000ffffffffffffffffbce6faada7179e84f3b9cac2fc632551, 0x1)
-    print(cmd_set_curve(group))
+    pass
+    # model = ShortWeierstrassModel()
+    # coords = model.coordinates["projective"]
+    # p = 0xffffffff00000001000000000000000000000000ffffffffffffffffffffffff
+    # curve = EllipticCurve(model, coords,
+    #                       p,
+    #                       {"a": 0xffffffff00000001000000000000000000000000fffffffffffffffffffffffc,
+    #                        "b": 0x5ac635d8aa3a93e7b3ebbd55769886bc651d06b0cc53b0f63bce3c3e27d2604b})
+    # affine_g = Point(AffineCoordinateModel(model),
+    #                  x=Mod(0x6b17d1f2e12c4247f8bce6e563a440f277037d812deb33a0f4a13945d898c296, p),
+    #                  y=Mod(0x4fe342e2fe1a7f9b8ee7eb4a7c0f9e162bce33576b315ececbb6406837bf51f5, p))
+    # g = Point.from_affine(coords, affine_g)
+    # neutral = Point(coords, X=Mod(0, p), Y=Mod(1, p), Z=Mod(0, p))
+    # group = AbelianGroup(curve, g, neutral,
+    #                      0xffffffff00000000ffffffffffffffffbce6faada7179e84f3b9cac2fc632551, 0x1)
+    #
+    # print(cmd_set_curve(group))
+    # mul = LTRMultiplier(coords.formulas["add-1998-cmo"], coords.formulas["dbl-1998-cmo"], coords.formulas["z"])
+    # mul.init(group, g)
+    # res = mul.multiply(0x2EF035DF6D0634C7422161D08BCC794B5312E042DDB32B0135A4DE6E6345A555)
+    # rx = Mod(0x77E3FF34C12571970845CBEB1BE0A79E3ECEE187510C2B8894BA800F8164C954, p)
+    # ry = Mod(0x408A6A05607F9ACA97BB9A34EA643B107AADE0C9BB5EDB930EADE3009666B9D1, p)
+    # rz = Mod(0xC66ECD687C335D63A7030434CA70351191BAFF1C206332EFEA39FA3003E91646, p)
+    #
+    # ox = Mod(0x3B1E7733E3250C97EB9D00AE0394F0768902DD337FEAAF7C4F6B9588462920DD, p)
+    # oy = Mod(0xBA718497596C964E77F9666506505B1E730EE69D254E85AD44727DFFB2C7063E, p)
+    # oz = Mod(0x0000000000000000000000000000000000000000000000000000000000000001, p)
+    #
+    # pt = Point(coords, X=rx, Y=ry, Z=rz)
+    # ot = Point(coords, X=ox, Y=oy, Z=oz)
+    # print(ot, res)
+    # print(pt.equals(res))
+    # print(ot.equals(res))
+    # print(ot == res)
 
 
 if __name__ == "__main__":
diff --git a/pyecsca/codegen/mult/mult.h b/pyecsca/codegen/mult.h
index d2817b7..b46684f 100644
--- a/pyecsca/codegen/mult/mult.h
+++ b/pyecsca/codegen/mult.h
@@ -3,9 +3,6 @@
 
 #include "defs.h"
 
-#define MULT_NONE			  0
-#define MULT_DOUBLE_AND_ADD   1
-
 void scalar_mult(bn_t *scalar, point_t *point, curve_t *curve, point_t *out);
 
 #endif //MULT_H_
\ No newline at end of file
diff --git a/pyecsca/codegen/mult/double_and_add.c b/pyecsca/codegen/mult/double_and_add.c
deleted file mode 100644
index 1cead25..0000000
--- a/pyecsca/codegen/mult/double_and_add.c
+++ /dev/null
@@ -1,19 +0,0 @@
-#include "mult.h"
-#include "point.h"
-
-void scalar_mult(bn_t *scalar, point_t *point, curve_t *curve, point_t *out) {
-	point_t *q = point_copy(point);
-	point_t *r = point_copy(curve->neutral);
-
-    int nbits = bn_bit_length(scalar);
-    for (int i = nbits; i >= 0; i--) {
-        point_dbl(r, curve, r);
-        if (bn_get_bit(scalar, i) == 1) {
-            point_add(q, r, curve, r);
-        }
-    }
-    //point_scl(r, curve, r);
-    point_set(r, out);
-    point_free(q);
-    point_free(r);
-}
\ No newline at end of file
diff --git a/pyecsca/codegen/mult/mult.c b/pyecsca/codegen/mult/mult.c
deleted file mode 100644
index d0f60ae..0000000
--- a/pyecsca/codegen/mult/mult.c
+++ /dev/null
@@ -1,6 +0,0 @@
-
-#include "mult.h"
-
-#if MULT == MULT_DOUBLE_AND_ADD
-#include "double_and_add.c"
-#endif
\ No newline at end of file
diff --git a/pyecsca/codegen/point.h b/pyecsca/codegen/point.h
index 61e7a7d..1b7e8e1 100644
--- a/pyecsca/codegen/point.h
+++ b/pyecsca/codegen/point.h
@@ -11,6 +11,8 @@ void point_set(const point_t *from, point_t *out);
 
 void point_free(point_t *point);
 
+bool point_equals(const point_t *one, const point_t *other);
+
 void point_to_affine(point_t *point, curve_t *curve, bn_t *out_x, bn_t *out_y);
 
 void point_from_affine(bn_t *x, bn_t *y, curve_t *curve, point_t *out);
diff --git a/pyecsca/codegen/templates/Makefile b/pyecsca/codegen/templates/Makefile
index 27402e4..9526640 100644
--- a/pyecsca/codegen/templates/Makefile
+++ b/pyecsca/codegen/templates/Makefile
@@ -1,14 +1,14 @@
 TARGET = pyecsca-codegen
 
-SRC += main.c bn/bn.c asn1/asn1.c hash/hash.c mult/mult.c prng/prng.c gen/point.c gen/curve.c
+SRC += main.c bn/bn.c asn1/asn1.c hash/hash.c prng/prng.c gen/point.c gen/curve.c gen/mult.c
 
 PLATFORM = {{ platform }}
 
-CDEFS += -DHASH={{ hash_type }} -DMULT={{ mult_algo }} -DMOD_RAND={{ mod_rand }}
+CDEFS += -DHASH={{ hash_type }} -DMOD_RAND={{ mod_rand }}
 
-MKDIR_LIST += hash prng mult asn1 bn gen
+MKDIR_LIST += hash prng asn1 bn gen
 
-EXTRAINCDIRS += hash prng mult asn1 bn gen tommath
+EXTRAINCDIRS += hash prng asn1 bn gen tommath
 
 LDFLAGS += tommath/libtommath-{{ platform }}.a
 
diff --git a/pyecsca/codegen/templates/formula_add.c b/pyecsca/codegen/templates/formula_add.c
index 971dd0a..05fd541 100644
--- a/pyecsca/codegen/templates/formula_add.c
+++ b/pyecsca/codegen/templates/formula_add.c
@@ -1,3 +1,13 @@
 void point_add(const point_t *one, const point_t *other, const curve_t *curve, point_t *out_one) {
+	{%- if short_circuit %}
+		if (point_equals(one, curve->neutral)) {
+			point_set(other, out_one);
+			return;
+		}
+		if (point_equals(other, curve->neutral)) {
+			point_set(one, out_one);
+			return;
+		}
+	{%- endif %}
 	{%- include "ops.c" %}
 }
\ No newline at end of file
diff --git a/pyecsca/codegen/templates/formula_dadd.c b/pyecsca/codegen/templates/formula_dadd.c
index 5d3bd05..7f4d2e5 100644
--- a/pyecsca/codegen/templates/formula_dadd.c
+++ b/pyecsca/codegen/templates/formula_dadd.c
@@ -1,3 +1,4 @@
 void point_dadd(const point_t *one, const point_t *other, const point_t *diff, const curve_t *curve, point_t *out_one) {
+	// TODO: short-circuits
 	{%- include "ops.c" %}
 }
diff --git a/pyecsca/codegen/templates/formula_dbl.c b/pyecsca/codegen/templates/formula_dbl.c
index 89543a1..bf91efb 100644
--- a/pyecsca/codegen/templates/formula_dbl.c
+++ b/pyecsca/codegen/templates/formula_dbl.c
@@ -1,3 +1,9 @@
 void point_dbl(const point_t *one, const curve_t *curve, point_t *out_one) {
+	{%- if short_circuit %}
+		if (point_equals(one, curve->neutral)) {
+			point_set(one, out_one);
+			return;
+		}
+	{%- endif %}
 	{%- include "ops.c" %}
 }
\ No newline at end of file
diff --git a/pyecsca/codegen/templates/formula_ladd.c b/pyecsca/codegen/templates/formula_ladd.c
index 73f2325..81bcf73 100644
--- a/pyecsca/codegen/templates/formula_ladd.c
+++ b/pyecsca/codegen/templates/formula_ladd.c
@@ -1,3 +1,4 @@
 void point_ladd(const point_t *one, const point_t *other, const point_t *diff, const curve_t *curve, point_t *out_one, point_t *out_other) {
+	// TODO: short-circuits
 	{%- include "ops.c" %}
 }
\ No newline at end of file
diff --git a/pyecsca/codegen/templates/formula_neg.c b/pyecsca/codegen/templates/formula_neg.c
index 156053c..fc790de 100644
--- a/pyecsca/codegen/templates/formula_neg.c
+++ b/pyecsca/codegen/templates/formula_neg.c
@@ -1,3 +1,9 @@
 void point_neg(const point_t *one, const curve_t *curve, point_t *out_one) {
+	{%- if short_circuit %}
+		if (point_equals(one, curve->neutral)) {
+			point_set(one, out_one);
+			return;
+		}
+	{%- endif %}
 	{%- include "ops.c" %}
 }
\ No newline at end of file
diff --git a/pyecsca/codegen/templates/formula_scl.c b/pyecsca/codegen/templates/formula_scl.c
index a9595eb..b5db4d7 100644
--- a/pyecsca/codegen/templates/formula_scl.c
+++ b/pyecsca/codegen/templates/formula_scl.c
@@ -1,3 +1,9 @@
 void point_scl(const point_t *one, const curve_t *curve, point_t *out_one) {
+	{%- if short_circuit %}
+		if (point_equals(one, curve->neutral)) {
+			point_set(one, out_one);
+			return;
+		}
+	{%- endif %}
 	{%- include "ops.c" %}
 }
\ No newline at end of file
diff --git a/pyecsca/codegen/templates/formula_tpl.c b/pyecsca/codegen/templates/formula_tpl.c
index debfd37..9b5df25 100644
--- a/pyecsca/codegen/templates/formula_tpl.c
+++ b/pyecsca/codegen/templates/formula_tpl.c
@@ -1,3 +1,9 @@
 void point_tpl(const point_t *one, const curve_t *curve, point_t *out_one) {
+	{%- if short_circuit %}
+		if (point_equals(one, curve->neutral)) {
+			point_set(one, out_one);
+			return;
+		}
+	{%- endif %}
 	{%- include "ops.c" %}
 }
\ No newline at end of file
diff --git a/pyecsca/codegen/templates/main.c b/pyecsca/codegen/templates/main.c
index 14a4e58..1c9491b 100644
--- a/pyecsca/codegen/templates/main.c
+++ b/pyecsca/codegen/templates/main.c
@@ -2,10 +2,10 @@
 #include "simpleserial/simpleserial.h"
 #include "asn1/asn1.h"
 #include "hash/hash.h"
-#include "mult/mult.h"
 #include "bn/bn.h"
 #include "prng/prng.h"
 #include "gen/defs.h"
+#include "mult.h"
 #include "point.h"
 #include "curve.h"
 #include "fat.h"
@@ -97,7 +97,7 @@ static uint8_t cmd_generate(uint8_t *data, uint16_t len) {
 	{%- for variable in curve_variables %}
 	bn_to_binpad(&pubkey->{{ variable }}, pub + coord_size * {{ loop.index0 }}, coord_size);
 	{%- endfor %}
-	simpleserial_put('w', coord_size * {{ curve_parameters | length }}, pub);
+	simpleserial_put('w', coord_size * {{ curve_variables | length }}, pub);
 	return 0;
 }
 
@@ -138,7 +138,7 @@ static void parse_scalar_mult(const char *path, const uint8_t *data, size_t len,
 }
 
 static uint8_t cmd_scalar_mult(uint8_t *data, uint16_t len) {
-	// perform base point scalar mult with supplied scalar, return affine point.
+	// perform base point scalar mult with supplied scalar.
 	bn_t scalar; bn_init(&scalar);
 	parse_data(data, len, "", parse_scalar_mult, (void *) &scalar);
 	size_t coord_size = bn_to_bin_size(&curve->p);
@@ -151,7 +151,7 @@ static uint8_t cmd_scalar_mult(uint8_t *data, uint16_t len) {
 	{%- for variable in curve_variables %}
 	bn_to_binpad(&result->{{ variable }}, res + coord_size * {{ loop.index0 }}, coord_size);
 	{%- endfor %}
-	simpleserial_put('w', coord_size * {{ curve_parameters | length }}, res);
+	simpleserial_put('w', coord_size * {{ curve_variables | length }}, res);
 	bn_clear(&scalar);
 	point_free(result);
 	return 0;
diff --git a/pyecsca/codegen/templates/mult.c b/pyecsca/codegen/templates/mult.c
new file mode 100644
index 0000000..61ec825
--- /dev/null
+++ b/pyecsca/codegen/templates/mult.c
@@ -0,0 +1,18 @@
+
+#include "mult.h"
+
+{%- if isinstance(scalarmult, LTRMultiplier) -%}
+{%- include "mult_ltr.c" %}
+{%- elif isinstance(scalarmult, RTLMultiplier) -%}
+{%- include "mult_rtl.c" %}
+{%- elif isinstance(scalarmult, CoronMultiplier) -%}
+{%- include "mult_coron.c" %}
+{%- elif isinstance(scalarmult, LadderMultiplier) -%}
+{%- include "mult_ldr.c" %}
+{%- elif isinstance(scalarmult, SimpleLadderMultiplier) -%}
+{%- include "mult_simple_ldr.c" %}
+{%- elif isinstance(scalarmult, DifferentialLadderMultiplier) -%}
+{%- include "mult_diff_ldr.c" %}
+{%- elif isinstance(scalarmult, BinaryNAFMultiplier) -%}
+{%- include "mult_bnaf.c" %}
+{%- endif -%}
diff --git a/pyecsca/codegen/templates/mult_bnaf.c b/pyecsca/codegen/templates/mult_bnaf.c
new file mode 100644
index 0000000..0b9c4de
--- /dev/null
+++ b/pyecsca/codegen/templates/mult_bnaf.c
@@ -0,0 +1,28 @@
+#include "mult.h"
+#include "point.h"
+
+void scalar_mult(bn_t *scalar, point_t *point, curve_t *curve, point_t *out) {
+	point_t *neg = point_new();
+	point_neg(point, curve, neg);
+	point_t *q = point_copy(&curve->neutral);
+
+	wnaf_t *naf = bn_bnaf(scalar);
+
+	for (size_t i = naf->length; i >= 0; i--) {
+		point_dbl(q, curve, q);
+		if (naf->data[i] == 1) {
+			point_add(q, point, curve, q);
+		} else if (naf->data[i] == -1) {
+			point_add(q, neg, curve, q);
+		}
+	}
+	free(naf->data);
+	free(naf);
+
+    {%- if "scl" in scalarmult.formulas %}
+    	point_scl(q, curve, q);
+    {%- endif %}
+    point_set(q, out);
+	point_free(neg);
+	point_free(q);
+}
\ No newline at end of file
diff --git a/pyecsca/codegen/templates/mult_coron.c b/pyecsca/codegen/templates/mult_coron.c
new file mode 100644
index 0000000..d8e8a0d
--- /dev/null
+++ b/pyecsca/codegen/templates/mult_coron.c
@@ -0,0 +1,22 @@
+#include "mult.h"
+#include "point.h"
+
+void scalar_mult(bn_t *scalar, point_t *point, curve_t *curve, point_t *out) {
+	point_t *p0 = point_copy(point);
+	point_t *p1 = point_new();
+
+	int nbits = bn_bit_length(scalar);
+	for (int i = nbits - 2; i >= 0; i--) {
+	  point_dbl(p0, curve, p0);
+	  point_add(p0, point, curve, p1);
+	  if (bn_get_bit(scalar, i) != 0) {
+		 point_set(p1, p0);
+	  }
+	}
+	{%- if "scl" in scalarmult.formulas %}
+		point_scl(p0, curve, p0);
+	{%- endif %}
+	point_set(p0, out);
+	point_free(p0);
+	point_free(p1);
+}
\ No newline at end of file
diff --git a/pyecsca/codegen/templates/mult_diff_ldr.c b/pyecsca/codegen/templates/mult_diff_ldr.c
new file mode 100644
index 0000000..2683116
--- /dev/null
+++ b/pyecsca/codegen/templates/mult_diff_ldr.c
@@ -0,0 +1,29 @@
+#include "mult.h"
+#include "point.h"
+
+void scalar_mult(bn_t *scalar, point_t *point, curve_t *curve, point_t *out) {
+	point_t *p0 = point_copy(&curve->neutral);
+	point_t *p1 = point_copy(point);
+	{%- if scalarmult.complete %}
+		int nbits = bn_bit_length(&curve->n) - 1;
+	{%- else %}
+		int nbits = bn_bit_length(scalar) - 1;
+	{%- endif %}
+
+	for (int i = nbits; i >= 0; i--) {
+		if (bn_get_bit(scalar, i) == 0) {
+			point_dadd(point, p0, p1, curve, p1);
+			point_dbl(p0, curve, p0);
+		} else {
+			point_dadd(point, p0, p1, curve, p0);
+			point_dbl(p1, curve, p1);
+		}
+	}
+
+	{%- if "scl" in scalarmult.formulas %}
+    	point_scl(p0, curve, p0);
+    {%- endif %}
+	point_set(p0, out);
+	point_free(p0);
+	point_free(p1);
+}
\ No newline at end of file
diff --git a/pyecsca/codegen/templates/mult_ldr.c b/pyecsca/codegen/templates/mult_ldr.c
new file mode 100644
index 0000000..cabab0b
--- /dev/null
+++ b/pyecsca/codegen/templates/mult_ldr.c
@@ -0,0 +1,30 @@
+#include "mult.h"
+#include "point.h"
+
+void scalar_mult(bn_t *scalar, point_t *point, curve_t *curve, point_t *out) {
+	{%- if scalarmult.complete %}
+		point_t *p0 = point_copy(&curve->neutral);
+		point_t *p1 = point_copy(point);
+		int nbits = bn_bit_length(&curve->n) - 1;
+	{%- else %}
+		point_t *p0 = point_copy(point);
+		point_t *p1 = point_new();
+		point_dbl(point, curve, p1);
+		int nbits = bn_bit_length(scalar) - 2;
+	{%- endif %}
+
+	for (int i = nbits; i >= 0; i--) {
+		if (bn_get_bit(scalar, i) == 0) {
+			point_ladd(p0, p1, curve, p0, p1);
+		} else {
+			point_ladd(p1, p0, curve, p1, p0);
+		}
+	}
+
+	{%- if "scl" in scalarmult.formulas %}
+    	point_scl(p0, curve, p0);
+    {%- endif %}
+	point_set(p0, out);
+	point_free(p0);
+	point_free(p1);
+}
\ No newline at end of file
diff --git a/pyecsca/codegen/templates/mult_ltr.c b/pyecsca/codegen/templates/mult_ltr.c
new file mode 100644
index 0000000..9b7d87e
--- /dev/null
+++ b/pyecsca/codegen/templates/mult_ltr.c
@@ -0,0 +1,40 @@
+#include "mult.h"
+#include "point.h"
+
+void scalar_mult(bn_t *scalar, point_t *point, curve_t *curve, point_t *out) {
+	{%- if scalarmult.complete %}
+		point_t *q = point_copy(point);
+		point_t *r = point_copy(curve->neutral);
+		int nbits = bn_bit_length(&curve->n) - 1;
+	{%- else %}
+		point_t *q = point_new();
+		point_dbl(point, curve, q);
+		point_t *r = point_copy(point);
+		int nbits = bn_bit_length(scalar) - 2;
+	{%- endif %}
+
+	{%- if scalarmult.always %}
+		point_t *dummy = point_new();
+	{%- endif %}
+
+    for (int i = nbits; i >= 0; i--) {
+        point_dbl(r, curve, r);
+        if (bn_get_bit(scalar, i) == 1) {
+            point_add(q, r, curve, r);
+        } else {
+        	{%- if scalarmult.always %}
+			point_add(q, r, curve, dummy);
+			{%- endif %}
+        }
+    }
+    {%- if "scl" in scalarmult.formulas %}
+    	point_scl(r, curve, r);
+    {%- endif %}
+
+    point_set(r, out);
+    point_free(q);
+    point_free(r);
+	{%- if scalarmult.always %}
+		point_free(dummy);
+	{%- endif %}
+}
\ No newline at end of file
diff --git a/pyecsca/codegen/templates/mult_rtl.c b/pyecsca/codegen/templates/mult_rtl.c
new file mode 100644
index 0000000..af437e0
--- /dev/null
+++ b/pyecsca/codegen/templates/mult_rtl.c
@@ -0,0 +1,37 @@
+#include "mult.h"
+#include "point.h"
+
+void scalar_mult(bn_t *scalar, point_t *point, curve_t *curve, point_t *out) {
+	point_t *q = point_copy(point);
+	point_t *r = point_copy(curve->neutral);
+
+	{%- if scalarmult.always %}
+		point_t *dummy = point_new();
+	{%- endif %}
+	bn_t copy;
+	bn_init(&copy);
+	bn_copy(scalar, &copy);
+
+    while (!bn_is_0(&copy)) {
+        if (bn_get_bit(&copy, i) == 1) {
+            point_add(q, r, curve, r);
+        } else {
+        	{%- if scalarmult.always %}
+			point_add(q, r, curve, dummy);
+			{%- endif %}
+        }
+        point_dbl(q, curve, q);
+        bn_rsh(&copy, 1, &copy);
+    }
+    {%- if "scl" in scalarmult.formulas %}
+    	point_scl(r, curve, r);
+    {%- endif %}
+
+    point_set(r, out);
+    point_free(q);
+    point_free(r);
+    bn_clear(&copy);
+	{%- if scalarmult.always %}
+		point_free(dummy);
+	{%- endif %}
+}
\ No newline at end of file
diff --git a/pyecsca/codegen/templates/mult_simple_ldr.c b/pyecsca/codegen/templates/mult_simple_ldr.c
new file mode 100644
index 0000000..438a44b
--- /dev/null
+++ b/pyecsca/codegen/templates/mult_simple_ldr.c
@@ -0,0 +1,29 @@
+#include "mult.h"
+#include "point.h"
+
+void scalar_mult(bn_t *scalar, point_t *point, curve_t *curve, point_t *out) {
+	point_t *p0 = point_copy(&curve->neutral);
+	point_t *p1 = point_copy(point);
+	{%- if scalarmult.complete %}
+		int nbits = bn_bit_length(&curve->n) - 1;
+	{%- else %}
+		int nbits = bn_bit_length(scalar) - 1;
+	{%- endif %}
+
+	for (int i = nbits; i >= 0; i--) {
+		if (bn_get_bit(scalar, i) == 1) {
+			point_add(p0, p1, curve, p1);
+			point_dbl(p0, curve, p0);
+		} else {
+			point_add(p0, p1, curve, p0);
+			point_dbl(p1, curve, p1);
+		}
+	}
+
+	{%- if "scl" in scalarmult.formulas %}
+    	point_scl(p0, curve, p0);
+    {%- endif %}
+	point_set(p0, out);
+	point_free(p0);
+	point_free(p1);
+}
\ No newline at end of file
diff --git a/pyecsca/codegen/templates/point.c b/pyecsca/codegen/templates/point.c
index 973fc42..689d598 100644
--- a/pyecsca/codegen/templates/point.c
+++ b/pyecsca/codegen/templates/point.c
@@ -29,6 +29,15 @@ void point_free(point_t *point) {
 	free(point);
 }
 
+bool point_equals(const point_t *one, const point_t *other) {
+	{%- for variable in variables %}
+	if (!bn_eq(&one->{{ variable }}, &other->{{ variable }})) {
+		return false;
+	}
+	{%- endfor %}
+	return true;
+}
+
 void point_to_affine(point_t *point, curve_t *curve, bn_t *out_x, bn_t *out_y) {
 	{%- include "ops.c" %}
 	{%- if "x" in allocations %}