Fix some scalar-mult algos, rename some outputs of formulas.

30 files changed, 131 insertions, 77 deletions

diff --git a/pyecsca/ec/efd/montgom/xz/diffadd/dadd-1987-m b/pyecsca/ec/efd/montgom/xz/diffadd/dadd-1987-m
index 31f105a..b0e2c79 100644
--- a/pyecsca/ec/efd/montgom/xz/diffadd/dadd-1987-m
+++ b/pyecsca/ec/efd/montgom/xz/diffadd/dadd-1987-m
@@ -1,3 +1,3 @@
 source 1987 Montgomery "Speeding the Pollard and elliptic curve methods of factorization", page 261, third display
-compute X5 = Z1(X2 X3-Z2 Z3)^2
-compute Z5 = X1(X2 Z3-Z2 X3)^2
+compute X4 = Z1(X2 X3-Z2 Z3)^2
+compute Z4 = X1(X2 Z3-Z2 X3)^2
diff --git a/pyecsca/ec/efd/montgom/xz/diffadd/dadd-1987-m-2 b/pyecsca/ec/efd/montgom/xz/diffadd/dadd-1987-m-2
index 814c28d..8d5ef46 100644
--- a/pyecsca/ec/efd/montgom/xz/diffadd/dadd-1987-m-2
+++ b/pyecsca/ec/efd/montgom/xz/diffadd/dadd-1987-m-2
@@ -1,3 +1,3 @@
 source 1987 Montgomery "Speeding the Pollard and elliptic curve methods of factorization", page 261, fifth display
-compute X5 = Z1((X3-Z3)(X2+Z2)+(X3+Z3)(X2-Z2))^2
-compute Z5 = X1((X3-Z3)(X2+Z2)-(X3+Z3)(X2-Z2))^2
+compute X4 = Z1((X3-Z3)(X2+Z2)+(X3+Z3)(X2-Z2))^2
+compute Z4 = X1((X3-Z3)(X2+Z2)-(X3+Z3)(X2-Z2))^2
diff --git a/pyecsca/ec/efd/montgom/xz/diffadd/dadd-1987-m-2.op3 b/pyecsca/ec/efd/montgom/xz/diffadd/dadd-1987-m-2.op3
index 2d31113..8cc9c2c 100644
--- a/pyecsca/ec/efd/montgom/xz/diffadd/dadd-1987-m-2.op3
+++ b/pyecsca/ec/efd/montgom/xz/diffadd/dadd-1987-m-2.op3
@@ -6,7 +6,7 @@ t4 = t2*t3
 t5 = t0*t1
 t6 = t5+t4
 t7 = t6^2
-X5 = Z1*t7
+X4 = Z1*t7
 t8 = X3-Z3
 t9 = X2+Z2
 t10 = X3+Z3
@@ -15,4 +15,4 @@ t12 = t10*t11
 t13 = t8*t9
 t14 = t13-t12
 t15 = t14^2
-Z5 = X1*t15
+Z4 = X1*t15
diff --git a/pyecsca/ec/efd/montgom/xz/diffadd/dadd-1987-m-3 b/pyecsca/ec/efd/montgom/xz/diffadd/dadd-1987-m-3
index 757dff5..adb5e88 100644
--- a/pyecsca/ec/efd/montgom/xz/diffadd/dadd-1987-m-3
+++ b/pyecsca/ec/efd/montgom/xz/diffadd/dadd-1987-m-3
@@ -5,5 +5,5 @@ compute C = X3+Z3
 compute D = X3-Z3
 compute DA = D A
 compute CB = C B
-compute X5 = Z1(DA+CB)^2
-compute Z5 = X1(DA-CB)^2
+compute X4 = Z1(DA+CB)^2
+compute Z4 = X1(DA-CB)^2
diff --git a/pyecsca/ec/efd/montgom/xz/diffadd/dadd-1987-m-3.op3 b/pyecsca/ec/efd/montgom/xz/diffadd/dadd-1987-m-3.op3
index 09bb8ef..011bb85 100644
--- a/pyecsca/ec/efd/montgom/xz/diffadd/dadd-1987-m-3.op3
+++ b/pyecsca/ec/efd/montgom/xz/diffadd/dadd-1987-m-3.op3
@@ -6,7 +6,7 @@ DA = D*A
 CB = C*B
 t0 = DA+CB
 t1 = t0^2
-X5 = Z1*t1
+X4 = Z1*t1
 t2 = DA-CB
 t3 = t2^2
-Z5 = X1*t3
+Z4 = X1*t3
diff --git a/pyecsca/ec/efd/montgom/xz/diffadd/dadd-1987-m.op3 b/pyecsca/ec/efd/montgom/xz/diffadd/dadd-1987-m.op3
index aa88ba8..5e55d59 100644
--- a/pyecsca/ec/efd/montgom/xz/diffadd/dadd-1987-m.op3
+++ b/pyecsca/ec/efd/montgom/xz/diffadd/dadd-1987-m.op3
@@ -2,9 +2,9 @@ t0 = Z2*Z3
 t1 = X2*X3
 t2 = t1-t0
 t3 = t2^2
-X5 = Z1*t3
+X4 = Z1*t3
 t4 = Z2*X3
 t5 = X2*Z3
 t6 = t5-t4
 t7 = t6^2
-Z5 = X1*t7
+Z4 = X1*t7
diff --git a/pyecsca/ec/efd/montgom/xz/diffadd/mdadd-1987-m b/pyecsca/ec/efd/montgom/xz/diffadd/mdadd-1987-m
index 5323f23..4ae695c 100644
--- a/pyecsca/ec/efd/montgom/xz/diffadd/mdadd-1987-m
+++ b/pyecsca/ec/efd/montgom/xz/diffadd/mdadd-1987-m
@@ -6,5 +6,5 @@ compute C = X3+Z3
 compute D = X3-Z3
 compute DA = D A
 compute CB = C B
-compute X5 = (DA+CB)^2
-compute Z5 = X1(DA-CB)^2
+compute X4 = (DA+CB)^2
+compute Z4 = X1(DA-CB)^2
diff --git a/pyecsca/ec/efd/montgom/xz/diffadd/mdadd-1987-m.op3 b/pyecsca/ec/efd/montgom/xz/diffadd/mdadd-1987-m.op3
index c734dfa..3b321d0 100644
--- a/pyecsca/ec/efd/montgom/xz/diffadd/mdadd-1987-m.op3
+++ b/pyecsca/ec/efd/montgom/xz/diffadd/mdadd-1987-m.op3
@@ -5,7 +5,7 @@ D = X3-Z3
 DA = D*A
 CB = C*B
 t0 = DA+CB
-X5 = t0^2
+X4 = t0^2
 t1 = DA-CB
 t2 = t1^2
-Z5 = X1*t2
+Z4 = X1*t2
diff --git a/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it b/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it
index 3a54cef..047a1f9 100644
--- a/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it
+++ b/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it
@@ -1,3 +1,3 @@
 source 2002 Izu--Takagi "A fast parallel elliptic curve multiplication resistant against side channel attacks", formula (8)
-compute X5 = Z1 ((X2 X3 - a Z2 Z3)^2 - 4 b Z2 Z3(X2 Z3 + X3 Z2))
-compute Z5 = X1 (X2 Z3 - X3 Z2)^2
+compute X4 = Z1 ((X2 X3 - a Z2 Z3)^2 - 4 b Z2 Z3(X2 Z3 + X3 Z2))
+compute Z4 = X1 (X2 Z3 - X3 Z2)^2
diff --git a/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it-2 b/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it-2
index 2280ca2..ca43635 100644
--- a/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it-2
+++ b/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it-2
@@ -1,5 +1,5 @@
 source 2002 Izu--Takagi "A fast parallel elliptic curve multiplication resistant against side channel attacks", formula (9)
 compute R = 2(X2 Z3 + X3 Z2)(X2 X3 + a Z2 Z3) + 4 b Z2^2 Z3^2
 compute S = (X2 Z3 - X3 Z2)^2
-compute X5 = R Z1 - S X1
-compute Z5 = S Z1
+compute X4 = R Z1 - S X1
+compute Z4 = S Z1
diff --git a/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it-2.op3 b/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it-2.op3
index f30ee85..3052fbc 100644
--- a/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it-2.op3
+++ b/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it-2.op3
@@ -19,5 +19,5 @@ t16 = t15-t14
 S = t16^2
 t17 = S*X1
 t18 = R*Z1
-X5 = t18-t17
-Z5 = S*Z1
+X4 = t18-t17
+Z4 = S*Z1
diff --git a/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it-3 b/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it-3
index 60a5382..f0fb438 100644
--- a/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it-3
+++ b/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it-3
@@ -11,7 +11,7 @@ compute T9 = 4 T8
 compute T10 = T3 + T4
 compute T11 = T9 T10
 compute T12 = T7 - T11
-compute X5 = Z1 T12
+compute X4 = Z1 T12
 compute T13 = T3 - T4
 compute T14 = T13^2
-compute Z5 = X1 T14
+compute Z4 = X1 T14
diff --git a/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it-3.op3 b/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it-3.op3
index cbfa800..0dc1e5e 100644
--- a/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it-3.op3
+++ b/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it-3.op3
@@ -10,7 +10,7 @@ T9 = 4*T8
 T10 = T3+T4
 T11 = T9*T10
 T12 = T7-T11
-X5 = Z1*T12
+X4 = Z1*T12
 T13 = T3-T4
 T14 = T13^2
-Z5 = X1*T14
+Z4 = X1*T14
diff --git a/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it-4 b/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it-4
index 978e636..9197ae8 100644
--- a/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it-4
+++ b/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it-4
@@ -16,5 +16,5 @@ compute T14 = T3 - T4
 compute T15 = T14^2
 compute T16 = Z1 T13
 compute T17 = X1 T15
-compute X5 = T16 - T17
-compute Z5 = Z1 T15
+compute X4 = T16 - T17
+compute Z4 = Z1 T15
diff --git a/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it-4.op3 b/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it-4.op3
index 33dda8c..bcd508e 100644
--- a/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it-4.op3
+++ b/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it-4.op3
@@ -15,5 +15,5 @@ T14 = T3-T4
 T15 = T14^2
 T16 = Z1*T13
 T17 = X1*T15
-X5 = T16-T17
-Z5 = Z1*T15
+X4 = T16-T17
+Z4 = Z1*T15
diff --git a/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it.op3 b/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it.op3
index 24f68c8..3e3f735 100644
--- a/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it.op3
+++ b/pyecsca/ec/efd/shortw/xz/diffadd/dadd-2002-it.op3
@@ -11,9 +11,9 @@ t9 = Z2*t8
 t10 = b*t9
 t11 = 4*t10
 t12 = t7-t11
-X5 = Z1*t12
+X4 = Z1*t12
 t13 = X3*Z2
 t14 = X2*Z3
 t15 = t14-t13
 t16 = t15^2
-Z5 = X1*t16
+Z4 = X1*t16
diff --git a/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-bj b/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-bj
index 6b19191..9539822 100644
--- a/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-bj
+++ b/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-bj
@@ -1,4 +1,4 @@
 source 2002 Brier--Joye "Weierstrass elliptic curves and side-channel attacks", formula (9) accompanied by note "7 multiplications plus 3 multiplications by a constant"
 assume Z1 = 1
-compute X5 = (X2 X3 - a Z2 Z3)^2 - 4 b Z2 Z3 (X2 Z3 + X3 Z2) 
-compute Z5 = X1 (X2 Z3 - X3 Z2)^2
+compute X4 = (X2 X3 - a Z2 Z3)^2 - 4 b Z2 Z3 (X2 Z3 + X3 Z2)
+compute Z4 = X1 (X2 Z3 - X3 Z2)^2
diff --git a/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-bj-2 b/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-bj-2
index fd943e2..33c0571 100644
--- a/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-bj-2
+++ b/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-bj-2
@@ -6,5 +6,5 @@ compute A = X2 X3
 compute B = Z2 Z3
 compute C = X2 Z3
 compute D = Z2 X3
-compute X5 = (A - a B)^2 - b4 B (C + D)
-compute Z5 = X1 (C - D)^2
+compute X4 = (A - a B)^2 - b4 B (C + D)
+compute Z4 = X1 (C - D)^2
diff --git a/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-bj-2.op3 b/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-bj-2.op3
index 10ec0ff..e8dcdf8 100644
--- a/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-bj-2.op3
+++ b/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-bj-2.op3
@@ -8,7 +8,7 @@ t2 = C+D
 t3 = t1^2
 t4 = B*t2
 t5 = b4*t4
-X5 = t3-t5
+X4 = t3-t5
 t6 = C-D
 t7 = t6^2
-Z5 = X1*t7
+Z4 = X1*t7
diff --git a/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-bj.op3 b/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-bj.op3
index c8e410d..496f313 100644
--- a/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-bj.op3
+++ b/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-bj.op3
@@ -10,9 +10,9 @@ t8 = Z3*t6
 t9 = Z2*t8
 t10 = b*t9
 t11 = 4*t10
-X5 = t7-t11
+X4 = t7-t11
 t12 = X3*Z2
 t13 = X2*Z3
 t14 = t13-t12
 t15 = t14^2
-Z5 = X1*t15
+Z4 = X1*t15
diff --git a/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it b/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it
index fe58ab2..4cdcc9f 100644
--- a/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it
+++ b/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it
@@ -1,4 +1,4 @@
 source 2002 Izu--Takagi "A fast parallel elliptic curve multiplication resistant against side channel attacks", formula (8), plus assumption Z1 = 1
 assume Z1 = 1
-compute X5 = (X2 X3 - a Z2 Z3)^2 - 4 b Z2 Z3(X2 Z3 + X3 Z2)
-compute Z5 = X1 (X2 Z3 - X3 Z2)^2
+compute X4 = (X2 X3 - a Z2 Z3)^2 - 4 b Z2 Z3(X2 Z3 + X3 Z2)
+compute Z4 = X1 (X2 Z3 - X3 Z2)^2
diff --git a/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it-2 b/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it-2
index a4ca509..6780268 100644
--- a/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it-2
+++ b/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it-2
@@ -2,5 +2,5 @@ source 2002 Izu--Takagi "A fast parallel elliptic curve multiplication resistant
 assume Z1 = 1
 compute R = 2(X2 Z3 + X3 Z2)(X2 X3 + a Z2 Z3) + 4 b Z2^2 Z3^2
 compute S = (X2 Z3 - X3 Z2)^2
-compute X5 = R - S X1
-compute Z5 = S
+compute X4 = R - S X1
+compute Z4 = S
diff --git a/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it-2.op3 b/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it-2.op3
index 8f2d95c..e3e1e4f 100644
--- a/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it-2.op3
+++ b/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it-2.op3
@@ -18,5 +18,5 @@ t15 = X2*Z3
 t16 = t15-t14
 S = t16^2
 t17 = S*X1
-X5 = R-t17
-Z5 = S
+X4 = R-t17
+Z4 = S
diff --git a/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it-3 b/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it-3
index 21daab3..3ed7f08 100644
--- a/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it-3
+++ b/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it-3
@@ -12,7 +12,7 @@ compute T9 = 4 T8
 compute T10 = T3 + T4
 compute T11 = T9 T10
 compute T12 = T7 - T11
-compute X5 = T12
+compute X4 = T12
 compute T13 = T3 - T4
 compute T14 = T13^2
-compute Z5 = X1 T14
+compute Z4 = X1 T14
diff --git a/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it-3.op3 b/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it-3.op3
index 01cffd0..f28599f 100644
--- a/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it-3.op3
+++ b/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it-3.op3
@@ -10,7 +10,7 @@ T9 = 4*T8
 T10 = T3+T4
 T11 = T9*T10
 T12 = T7-T11
-X5 = T12
+X4 = T12
 T13 = T3-T4
 T14 = T13^2
-Z5 = X1*T14
+Z4 = X1*T14
diff --git a/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it-4 b/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it-4
index 839c746..add2e0f 100644
--- a/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it-4
+++ b/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it-4
@@ -17,5 +17,5 @@ compute T14 = T3 - T4
 compute T15 = T14^2
 compute T16 = T13
 compute T17 = X1 T15
-compute X5 = T16 - T17
-compute Z5 = T15
+compute X4 = T16 - T17
+compute Z4 = T15
diff --git a/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it-4.op3 b/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it-4.op3
index 5363020..ef831fc 100644
--- a/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it-4.op3
+++ b/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it-4.op3
@@ -15,5 +15,5 @@ T14 = T3-T4
 T15 = T14^2
 T16 = T13
 T17 = X1*T15
-X5 = T16-T17
-Z5 = T15
+X4 = T16-T17
+Z4 = T15
diff --git a/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it.op3 b/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it.op3
index c8e410d..496f313 100644
--- a/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it.op3
+++ b/pyecsca/ec/efd/shortw/xz/diffadd/mdadd-2002-it.op3
@@ -10,9 +10,9 @@ t8 = Z3*t6
 t9 = Z2*t8
 t10 = b*t9
 t11 = 4*t10
-X5 = t7-t11
+X4 = t7-t11
 t12 = X3*Z2
 t13 = X2*Z3
 t14 = t13-t12
 t15 = t14^2
-Z5 = X1*t15
+Z4 = X1*t15
diff --git a/pyecsca/ec/mult.py b/pyecsca/ec/mult.py
index aba2739..4ccbabf 100644
--- a/pyecsca/ec/mult.py
+++ b/pyecsca/ec/mult.py
@@ -1,12 +1,12 @@
 from copy import copy
 from public import public
-from typing import Mapping, Tuple, Optional, MutableMapping
+from typing import Mapping, Tuple, Optional, MutableMapping, Union
 
 from pyecsca.ec.naf import naf, wnaf
 from .context import Context
 from .curve import EllipticCurve
 from .formula import (Formula, AdditionFormula, DoublingFormula, ScalingFormula, LadderFormula,
-                      NegationFormula)
+                      NegationFormula, DifferentialAdditionFormula)
 from .point import Point
 
 
@@ -54,6 +54,16 @@ class ScalarMultiplier(object):
         return self.context.execute(self.formulas["ladd"], start, to_dbl, to_add,
                                     **self.curve.parameters)
 
+    def _dadd(self, start: Point, one: Point, other: Point) -> Point:
+        if "dadd" not in self.formulas:
+            raise NotImplementedError
+        if one == self.curve.neutral:
+            return copy(other)
+        if other == self.curve.neutral:
+            return copy(one)
+        return self.context.execute(self.formulas["dadd"], start, one, other,
+                                    **self.curve.parameters)[0]
+
     def _neg(self, point: Point) -> Point:
         if "neg" not in self.formulas:
             raise NotImplementedError
@@ -91,9 +101,11 @@ class LTRMultiplier(ScalarMultiplier):
         self.always = always
 
     def multiply(self, scalar: int, point: Optional[Point] = None) -> Point:
+        if scalar == 0:
+            return copy(self.curve.neutral)
         q = self._init_multiply(point)
         r = copy(self.curve.neutral)
-        for i in range(scalar.bit_length(), -1, -1):
+        for i in range(scalar.bit_length() - 1, -1, -1):
             r = self._dbl(r)
             if scalar & (1 << i) != 0:
                 r = self._add(r, q)
@@ -120,6 +132,8 @@ class RTLMultiplier(ScalarMultiplier):
         self.always = always
 
     def multiply(self, scalar: int, point: Optional[Point] = None) -> Point:
+        if scalar == 0:
+            return copy(self.curve.neutral)
         q = self._init_multiply(point)
         r = copy(self.curve.neutral)
         while scalar > 0:
@@ -171,11 +185,13 @@ class LadderMultiplier(ScalarMultiplier):
         super().__init__(curve, ctx, ladd=ladd, scl=scl)
 
     def multiply(self, scalar: int, point: Optional[Point] = None) -> Point:
+        if scalar == 0:
+            return copy(self.curve.neutral)
         q = self._init_multiply(point)
-        p0 = copy(q)
-        p1 = self._ladd(self.curve.neutral, q, q)[1]
+        p0 = copy(self.curve.neutral)
+        p1 = copy(q)
         for i in range(scalar.bit_length() - 1, -1, -1):
-            if scalar & (1 << i) != 0:
+            if scalar & (1 << i) == 0:
                 p0, p1 = self._ladd(q, p0, p1)
             else:
                 p1, p0 = self._ladd(q, p1, p0)
@@ -189,22 +205,38 @@ class SimpleLadderMultiplier(ScalarMultiplier):
     """
     Montgomery ladder multiplier, using addition and doubling formulas.
     """
+    _differential: bool = False
 
-    def __init__(self, curve: EllipticCurve, add: AdditionFormula, dbl: DoublingFormula,
+    def __init__(self, curve: EllipticCurve,
+                 add: Union[AdditionFormula, DifferentialAdditionFormula], dbl: DoublingFormula,
                  scl: ScalingFormula = None, ctx: Context = None):
-        super().__init__(curve, ctx, add=add, dbl=dbl, scl=scl)
+        if isinstance(add, AdditionFormula):
+            super().__init__(curve, ctx, add=add, dbl=dbl, scl=scl)
+        elif isinstance(add, DifferentialAdditionFormula):
+            super().__init__(curve, ctx, dadd=add, dbl=dbl, scl=scl)
+            self._differential = True
+        else:
+            raise ValueError
 
     def multiply(self, scalar: int, point: Optional[Point] = None) -> Point:
+        if scalar == 0:
+            return copy(self.curve.neutral)
         q = self._init_multiply(point)
-        p0 = copy(q)
-        p1 = self._dbl(q)
-        for i in range(scalar.bit_length() - 2, -1, -1):
-            if scalar & (1 << i) != 0:
-                p0 = self._add(p0, p1)
-                p1 = self._dbl(p1)
-            else:
-                p1 = self._add(p0, p1)
+        p0 = copy(self.curve.neutral)
+        p1 = copy(q)
+        for i in range(scalar.bit_length() - 1, -1, -1):
+            if scalar & (1 << i) == 0:
+                if self._differential:
+                    p1 = self._dadd(q, p0, p1)
+                else:
+                    p1 = self._add(p0, p1)
                 p0 = self._dbl(p0)
+            else:
+                if self._differential:
+                    p0 = self._dadd(q, p0, p1)
+                else:
+                    p0 = self._add(p0, p1)
+                p1 = self._dbl(p1)
         if "scl" in self.formulas:
             p0 = self._scl(p0)
         return p0
@@ -226,6 +258,8 @@ class BinaryNAFMultiplier(ScalarMultiplier):
         self._point_neg = self._neg(point)
 
     def multiply(self, scalar: int, point: Optional[Point] = None) -> Point:
+        if scalar == 0:
+            return copy(self.curve.neutral)
         self._init_multiply(point)
         bnaf = naf(scalar)
         q = copy(self.curve.neutral)
@@ -246,23 +280,32 @@ class WindowNAFMultiplier(ScalarMultiplier):
     Window NAF (Non Adjacent Form) multiplier, left-to-right.
     """
     _points: MutableMapping[int, Point]
+    _points_neg: MutableMapping[int, Point]
+    _precompute_neg: bool = False
     _width: int
 
     def __init__(self, curve: EllipticCurve, add: AdditionFormula, dbl: DoublingFormula,
-                 neg: NegationFormula, width: int, scl: ScalingFormula = None, ctx: Context = None):
+                 neg: NegationFormula, width: int, scl: ScalingFormula = None, ctx: Context = None,
+                 precompute_negation: bool = False):
         super().__init__(curve, ctx, add=add, dbl=dbl, neg=neg, scl=scl)
         self._width = width
+        self._precompute_neg = precompute_negation
 
     def init(self, point: Point):
         self._point = point
         self._points = {}
+        self._points_neg = {}
         current_point = point
         double_point = self._dbl(point)
         for i in range(1, (self._width + 1) // 2 + 1):
             self._points[2 ** i - 1] = current_point
+            if self._precompute_neg:
+                self._points_neg[2 ** i - 1] = self._neg(current_point)
             current_point = self._add(current_point, double_point)
 
     def multiply(self, scalar: int, point: Optional[Point] = None):
+        if scalar == 0:
+            return copy(self.curve.neutral)
         self._init_multiply(point)
         naf = wnaf(scalar, self._width)
         q = copy(self.curve.neutral)
@@ -271,7 +314,10 @@ class WindowNAFMultiplier(ScalarMultiplier):
             if val > 0:
                 q = self._add(q, self._points[val])
             elif val < 0:
-                neg = self._neg(self._points[-val])
+                if self._precompute_neg:
+                    neg = self._points_neg[-val]
+                else:
+                    neg = self._neg(self._points[-val])
                 q = self._add(q, neg)
         if "scl" in self.formulas:
             q = self._scl(q)
diff --git a/pyecsca/ec/point.py b/pyecsca/ec/point.py
index 05f2503..0b9e9b3 100644
--- a/pyecsca/ec/point.py
+++ b/pyecsca/ec/point.py
@@ -1,5 +1,6 @@
 from typing import Mapping
 
+from pyecsca.ec.formula import ScalingFormula
 from .coordinates import CoordinateModel
 from .mod import Mod
 
@@ -18,8 +19,15 @@ class Point(object):
         # TODO: Somehow compare projective points. Via a map to an affinepoint?
         if type(other) is not Point:
             return False
-        return self.coordinate_model == other.coordinate_model and self.coords == other.coords
+        if  self.coordinate_model != other.coordinate_model:
+            return False
+        self_scaling = list(filter(lambda x: isinstance(x, ScalingFormula), self.coordinate_model.formulas.items()))
+        other_scaling = list(filter(lambda x: isinstance(x, ScalingFormula), other.coordinate_model.formulas.items()))
+        return self.coords == other.coords
 
-    def __repr__(self):
+    def __str__(self):
         args = ", ".join([f"{key}={val}" for key, val in self.coords.items()])
-        return f"Point([{args}] in {self.coordinate_model})"
+        return f"[{args}]"
+
+    def __repr__(self):
+        return f"Point([{str(self)}] in {self.coordinate_model})"