diff --git a/LICENSE-APACHE b/LICENSE-APACHE
index 91ba015..7774c4e 100644
--- a/LICENSE-APACHE
+++ b/LICENSE-APACHE
@@ -189,7 +189,7 @@ APPENDIX: How to apply the Apache License to your work.
 Copyright 2006-2009 Graydon Hoare
 Copyright 2009-2013 Mozilla Foundation
 Copyright 2018 Input Output HK
-Copyright 2018-2021 Typed IO
+Copyright 2018-2022 Typed IO
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
diff --git a/LICENSE-MIT b/LICENSE-MIT
index ec7d3ec..d944c2a 100644
--- a/LICENSE-MIT
+++ b/LICENSE-MIT
@@ -1,7 +1,7 @@
 Copyright (c) 2006-2009 Graydon Hoare
 Copyright (c) 2009-2013 Mozilla Foundation
 Copyright (c) 2018 Input Output HK
-Copyright (c) 2018-2021 Typed IO
+Copyright (c) 2018-2022 Typed IO
 
 Permission is hereby granted, free of charge, to any
 person obtaining a copy of this software and associated
diff --git a/src/curve25519/fe/fe32/mod.rs b/src/curve25519/fe/fe32/mod.rs
index fbbf36c..843f2ed 100644
--- a/src/curve25519/fe/fe32/mod.rs
+++ b/src/curve25519/fe/fe32/mod.rs
@@ -552,6 +552,19 @@ impl Fe {
             h5 as i32, h6 as i32, h7 as i32, h8 as i32, h9 as i32])
     }
 
+    pub(crate) fn negate_mut(&mut self) {
+        self.0[0] = -self.0[0];
+        self.0[1] = -self.0[1];
+        self.0[2] = -self.0[2];
+        self.0[3] = -self.0[3];
+        self.0[4] = -self.0[4];
+        self.0[5] = -self.0[5];
+        self.0[6] = -self.0[6];
+        self.0[7] = -self.0[7];
+        self.0[8] = -self.0[8];
+        self.0[9] = -self.0[9];
+    }
+
     /*
     h = f * f
     Can overlap h with f.
diff --git a/src/curve25519/fe/fe64/mod.rs b/src/curve25519/fe/fe64/mod.rs
index e56ece3..def15f8 100644
--- a/src/curve25519/fe/fe64/mod.rs
+++ b/src/curve25519/fe/fe64/mod.rs
@@ -10,6 +10,10 @@ use core::ops::{Add, Mul, Neg, Sub};
 
 pub mod precomp;
 
+// multiple of P
+const FOUR_P0: u64 = 0x1fffffffffffb4;
+const FOUR_P1234: u64 = 0x1ffffffffffffc;
+
 /// Field Element in \Z/(2^255-19)
 #[derive(Clone)]
 pub struct Fe(pub(crate) [u64; 5]);
@@ -90,10 +94,6 @@ impl Sub for &Fe {
 
     #[rustfmt::skip]
     fn sub(self, rhs: &Fe) -> Fe {
-        // multiple of P
-        const FOUR_P0: u64 = 0x1fffffffffffb4;
-        const FOUR_P1234: u64 = 0x1ffffffffffffc;
-
         let Fe([f0, f1, f2, f3, f4]) = *self;
         let Fe([g0, g1, g2, g3, g4]) = *rhs;
 
@@ -110,8 +110,17 @@ impl Sub for &Fe {
 impl Neg for &Fe {
     type Output = Fe;
 
+    #[rustfmt::skip]
     fn neg(self) -> Fe {
-        &Fe::ZERO - &self
+        let Fe([g0, g1, g2, g3, g4]) = *self;
+
+        let mut h0 = FOUR_P0    - g0    ; let c = h0 >> 51; h0 &= MASK;
+        let mut h1 = FOUR_P1234 - g1 + c; let c = h1 >> 51; h1 &= MASK;
+        let mut h2 = FOUR_P1234 - g2 + c; let c = h2 >> 51; h2 &= MASK;
+        let mut h3 = FOUR_P1234 - g3 + c; let c = h3 >> 51; h3 &= MASK;
+        let mut h4 = FOUR_P1234 - g4 + c; let c = h4 >> 51; h4 &= MASK;
+        h0 += c * 19;
+        Fe([h0, h1, h2, h3, h4])
     }
 }
 
@@ -279,6 +288,16 @@ impl Fe {
         Fe([r0, r1, r2, r3, r4])
     }
 
+    #[rustfmt::skip]
+    pub(crate) fn negate_mut(&mut self) {
+        self.0[0] = FOUR_P0    - self.0[0]    ; let c = self.0[0] >> 51; self.0[0] &= MASK;
+        self.0[1] = FOUR_P1234 - self.0[1] + c; let c = self.0[1] >> 51; self.0[1] &= MASK;
+        self.0[2] = FOUR_P1234 - self.0[2] + c; let c = self.0[2] >> 51; self.0[2] &= MASK;
+        self.0[3] = FOUR_P1234 - self.0[3] + c; let c = self.0[3] >> 51; self.0[3] &= MASK;
+        self.0[4] = FOUR_P1234 - self.0[4] + c; let c = self.0[4] >> 51; self.0[4] &= MASK;
+        self.0[0] += c * 19;
+    }
+
     /// Compute the square of the field element
     #[rustfmt::skip]
     pub fn square(&self) -> Fe {
diff --git a/src/curve25519/ge.rs b/src/curve25519/ge.rs
index 377bcd3..cae75d2 100644
--- a/src/curve25519/ge.rs
+++ b/src/curve25519/ge.rs
@@ -5,19 +5,28 @@ use super::fe::{precomp, Fe};
 use super::scalar::Scalar;
 use crate::constant_time::{Choice, CtEqual, CtZero};
 
+/// Curve Group Element (Point)
+///
+/// The group element is using the extended homogeneous coordinates
+/// using (x,y,z,t) which maps to coordinates using the following
+/// equations:
+///    X     = x/z
+///    Y     = y/z
+///    X * Y = t/z
 #[derive(Clone)]
-pub struct GeP2 {
+pub struct Ge {
     x: Fe,
     y: Fe,
     z: Fe,
+    t: Fe,
 }
 
+/// Curve Group element without t=X*Y
 #[derive(Clone)]
-pub struct GeP3 {
+pub struct GePartial {
     x: Fe,
     y: Fe,
     z: Fe,
-    t: Fe,
 }
 
 #[derive(Clone)]
@@ -35,6 +44,12 @@ pub struct GePrecomp {
     pub(crate) xy2d: Fe,
 }
 
+#[derive(Clone)]
+struct GeAffine {
+    x: Fe,
+    y: Fe,
+}
+
 #[derive(Clone)]
 pub struct GeCached {
     y_plus_x: Fe,
@@ -43,17 +58,62 @@ pub struct GeCached {
     t2d: Fe,
 }
 
+impl GeAffine {
+    pub fn to_bytes(&self) -> [u8; 32] {
+        let mut bs = self.y.to_bytes();
+        bs[31] ^= (if self.x.is_negative() { 1 } else { 0 }) << 7;
+        bs
+    }
+
+    pub fn from_bytes(s: &[u8; 32]) -> Option<Self> {
+        // See RFC8032 5.3.1 decoding process
+        //
+        // y (255 bits) | sign(x) (1 bit) = s
+        // let u = y^2 - 1
+        //     v = d * y^2 + 1
+        //     x = u * v^3 * (u * v^7)^((p-5)/8)
+
+        // recover y by clearing the highest bit (side effect of from_bytes)
+        let y = Fe::from_bytes(s);
+
+        // recover x
+        let y2 = y.square();
+        let u = &y2 - &Fe::ONE;
+        let v = &(&y2 * &Fe::D) + &Fe::ONE;
+        let v3 = &v.square() * &v;
+        let v7 = &v3.square() * &v;
+        let uv7 = &v7 * &u;
+
+        let mut x = &(&uv7.pow25523() * &v3) * &u;
+
+        let vxx = &x.square() * &v;
+        let check = &vxx - &u;
+        if check.is_nonzero() {
+            let check2 = &vxx + &u;
+            if check2.is_nonzero() {
+                return None;
+            }
+            x = &x * &Fe::SQRTM1;
+        }
+
+        if x.is_negative() == ((s[31] >> 7) != 0) {
+            x.negate_mut();
+        }
+        Some(Self { x, y })
+    }
+}
+
 impl GeP1P1 {
-    pub(crate) fn to_p2(&self) -> GeP2 {
-        GeP2 {
+    pub fn to_partial(&self) -> GePartial {
+        GePartial {
             x: &self.x * &self.t,
             y: &self.y * &self.z,
             z: &self.z * &self.t,
         }
     }
 
-    pub(crate) fn to_p3(&self) -> GeP3 {
-        GeP3 {
+    pub fn to_full(&self) -> Ge {
+        Ge {
             x: &self.x * &self.t,
             y: &self.y * &self.z,
             z: &self.z * &self.t,
@@ -62,8 +122,8 @@ impl GeP1P1 {
     }
 }
 
-impl GeP2 {
-    pub const ZERO: GeP2 = GeP2 {
+impl GePartial {
+    pub const ZERO: Self = Self {
         x: Fe::ZERO,
         y: Fe::ONE,
         z: Fe::ONE,
@@ -78,7 +138,7 @@ impl GeP2 {
         bs
     }
 
-    pub fn dbl(&self) -> GeP1P1 {
+    pub fn double_p1p1(&self) -> GeP1P1 {
         let xx = self.x.square();
         let yy = self.y.square();
         let b = self.z.square_and_double();
@@ -97,29 +157,49 @@ impl GeP2 {
         }
     }
 
-    /*
-    r = a * A + b * B
-    where a = a[0]+256*a[1]+...+256^31 a[31].
-    and b = b[0]+256*b[1]+...+256^31 b[31].
-    B is the Ed25519 base point (x,4/5) with x positive.
-    */
-    pub fn double_scalarmult_vartime(a_scalar: &Scalar, a_point: GeP3, b_scalar: &Scalar) -> GeP2 {
+    pub fn double(&self) -> Self {
+        self.double_p1p1().to_partial()
+    }
+
+    pub fn double_full(&self) -> Ge {
+        self.double_p1p1().to_full()
+    }
+
+    /// Calculate r = a * A + b * B
+    ///
+    /// ```ignore
+    /// double_scalarmult_vartime(a, A, b) = a * A + b * B
+    /// ```
+    ///
+    /// where
+    ///     a is a scalar
+    ///     A is an arbitrary point
+    ///     b is a scalar
+    ///     B the ED25519 base point (not a parameter to the function)
+    ///
+    /// Note that the
+    ///
+    pub fn double_scalarmult_vartime(
+        a_scalar: &Scalar,
+        a_point: Ge,
+        b_scalar: &Scalar,
+    ) -> GePartial {
         let aslide = a_scalar.slide();
         let bslide = b_scalar.slide();
 
         let a1 = a_point.to_cached();
-        let a2 = a_point.dbl().to_p3();
-        let a3 = (&a2 + &a1).to_p3().to_cached();
-        let a5 = (&a2 + &a3).to_p3().to_cached();
-        let a7 = (&a2 + &a5).to_p3().to_cached();
-        let a9 = (&a2 + &a7).to_p3().to_cached();
-        let a11 = (&a2 + &a9).to_p3().to_cached();
-        let a13 = (&a2 + &a11).to_p3().to_cached();
-        let a15 = (&a2 + &a13).to_p3().to_cached();
+        let a2 = a_point.double_p1p1().to_full();
+        let a3 = (&a2 + &a1).to_full().to_cached();
+        let a5 = (&a2 + &a3).to_full().to_cached();
+        let a7 = (&a2 + &a5).to_full().to_cached();
+        let a9 = (&a2 + &a7).to_full().to_cached();
+        let a11 = (&a2 + &a9).to_full().to_cached();
+        let a13 = (&a2 + &a11).to_full().to_cached();
+        let a15 = (&a2 + &a13).to_full().to_cached();
 
         let ai: [GeCached; 8] = [a1, a3, a5, a7, a9, a11, a13, a15];
 
-        let mut r = GeP2::ZERO;
+        let mut r = GePartial::ZERO;
 
         let mut i: usize = 255;
         loop {
@@ -133,20 +213,20 @@ impl GeP2 {
         }
 
         loop {
-            let mut t = r.dbl();
+            let mut t = r.double_p1p1();
             match aslide[i].cmp(&0) {
-                Ordering::Greater => t = &t.to_p3() + &ai[(aslide[i] / 2) as usize],
-                Ordering::Less => t = &t.to_p3() - &ai[(-aslide[i] / 2) as usize],
+                Ordering::Greater => t = &t.to_full() + &ai[(aslide[i] / 2) as usize],
+                Ordering::Less => t = &t.to_full() - &ai[(-aslide[i] / 2) as usize],
                 Ordering::Equal => {}
             }
 
             match bslide[i].cmp(&0) {
-                Ordering::Greater => t = &t.to_p3() + &precomp::BI[(bslide[i] / 2) as usize],
-                Ordering::Less => t = &t.to_p3() - &precomp::BI[(-bslide[i] / 2) as usize],
+                Ordering::Greater => t = &t.to_full() + &precomp::BI[(bslide[i] / 2) as usize],
+                Ordering::Less => t = &t.to_full() - &precomp::BI[(-bslide[i] / 2) as usize],
                 Ordering::Equal => {}
             }
 
-            r = t.to_p2();
+            r = t.to_partial();
 
             if i == 0 {
                 return r;
@@ -156,53 +236,51 @@ impl GeP2 {
     }
 }
 
-impl GeP3 {
-    pub fn from_bytes_negate_vartime(s: &[u8; 32]) -> Option<GeP3> {
-        // See RFC8032 5.3.1 decoding process
-        //
-        // y (255 bits) | sign(x) (1 bit) = s
-        // let u = y^2 - 1
-        //     v = d * y^2 + 1
-        //     x = u * v^3 * (u * v^7)^((p-5)/8)
-        let y = Fe::from_bytes(s);
-        let y2 = y.square();
-        let u = &y2 - &Fe::ONE;
-        let v = &(&y2 * &Fe::D) + &Fe::ONE;
-        let v3 = &v.square() * &v;
-        let v7 = &v3.square() * &v;
-        let uv7 = &v7 * &u;
-
-        let mut x = &(&uv7.pow25523() * &v3) * &u;
-
-        let vxx = &x.square() * &v;
-        let check = &vxx - &u;
-        if check.is_nonzero() {
-            let check2 = &vxx + &u;
-            if check2.is_nonzero() {
-                return None;
-            }
-            x = &x * &Fe::SQRTM1;
-        }
+impl Ge {
+    /// The Identity Element for the group, which represent (X=0, Y=1) and is (x=0, y=1, z=1, t=0*1)
+    pub const ZERO: Ge = Ge {
+        x: Fe::ZERO,
+        y: Fe::ONE,
+        z: Fe::ONE,
+        t: Fe::ZERO,
+    };
 
-        if x.is_negative() == ((s[31] >> 7) != 0) {
-            x = x.neg();
+    /// Create a group element from affine coordinate
+    #[inline]
+    fn from_affine(affine: GeAffine) -> Ge {
+        let t = &affine.x * &affine.y;
+        Ge {
+            x: affine.x,
+            y: affine.y,
+            z: Fe::ONE,
+            t: t,
         }
+    }
 
-        let t = &x * &y;
+    /// Flatten a group element on the affine plane (x,y)
+    #[inline]
+    fn to_affine(&self) -> GeAffine {
+        let recip = self.z.invert();
+        let x = &self.x * &recip;
+        let y = &self.y * &recip;
+        GeAffine { x, y }
+    }
 
-        Some(GeP3 {
-            x: x,
-            y: y,
-            z: Fe::ONE,
-            t: t,
-        })
+    /// Try to construct a group element (Point on the curve)
+    /// from its compressed byte representation (32 bytes little endian).
+    ///
+    /// The compressed bytes representation is the y coordinate (255 bits)
+    /// and the sign of the x coordinate (1 bit) as the highest bit.
+    pub fn from_bytes(s: &[u8; 32]) -> Option<Ge> {
+        GeAffine::from_bytes(s).map(Self::from_affine)
     }
 
-    pub fn to_p2(&self) -> GeP2 {
-        GeP2 {
-            x: self.x.clone(),
-            y: self.y.clone(),
-            z: self.z.clone(),
+    /// Drop the t coordinate to become a `GePartial`
+    pub fn to_partial(self) -> GePartial {
+        GePartial {
+            x: self.x,
+            y: self.y,
+            z: self.z,
         }
     }
 
@@ -215,35 +293,88 @@ impl GeP3 {
         }
     }
 
-    pub const ZERO: GeP3 = GeP3 {
-        x: Fe::ZERO,
-        y: Fe::ONE,
-        z: Fe::ONE,
-        t: Fe::ZERO,
-    };
+    pub fn double_p1p1(&self) -> GeP1P1 {
+        let xx = self.x.square();
+        let yy = self.y.square();
+        let b = self.z.square_and_double();
+        let a = &self.x + &self.y;
+        let aa = a.square();
+        let y3 = &yy + &xx;
+        let z3 = &yy - &xx;
+        let x3 = &aa - &y3;
+        let t3 = &b - &z3;
+
+        GeP1P1 {
+            x: x3,
+            y: y3,
+            z: z3,
+            t: t3,
+        }
+    }
 
-    pub fn dbl(&self) -> GeP1P1 {
-        self.to_p2().dbl()
+    /// Double the point
+    pub fn double(&self) -> Ge {
+        self.double_p1p1().to_full()
     }
 
+    pub fn double_partial(&self) -> GePartial {
+        self.double_p1p1().to_partial()
+    }
+
+    /// Transform a point into the compressed byte representation
+    ///
+    /// the compressed bytes representation is the Y coordinate
+    /// followed by the 1 bit sign of X coordinate
     pub fn to_bytes(&self) -> [u8; 32] {
-        let recip = self.z.invert();
-        let x = &self.x * &recip;
-        let y = &self.y * &recip;
-        let mut bs = y.to_bytes();
-        bs[31] ^= (if x.is_negative() { 1 } else { 0 }) << 7;
-        bs
+        self.to_affine().to_bytes()
     }
-}
 
-impl Add<GeCached> for GeP3 {
-    type Output = GeP1P1;
-    fn add(self, rhs: GeCached) -> GeP1P1 {
-        &self + &rhs
+    /// Compute r = a * B
+    ///
+    /// where
+    ///     a = a[0]+2^8*a[1]+...+2^248 a[31] a scalar number represented by 32-bytes in little endian format
+    ///         and a[31] is <= 0x80
+    ///     B the ED25519 base point (not a parameter to the function)
+    pub fn scalarmult_base(a: &Scalar) -> Ge {
+        let mut r: GeP1P1;
+        let mut t: GePrecomp;
+
+        /* each es[i] is between 0 and 0xf */
+        /* es[63] is between 0 and 7 */
+        let mut es = a.nibbles();
+
+        let mut carry: i8 = 0;
+        for esi in es[0..63].iter_mut() {
+            *esi += carry;
+            carry = *esi + 8;
+            carry >>= 4;
+            *esi -= carry << 4;
+        }
+        es[63] += carry;
+        /* each es[i] is between -8 and 8 */
+
+        let mut h = Ge::ZERO;
+        for j in 0..32 {
+            let i = j * 2 + 1;
+            t = GePrecomp::select(j, es[i]);
+            r = &h + &t;
+            h = r.to_full();
+        }
+
+        h = h.double_partial().double().double().double_full();
+
+        for j in 0..32 {
+            let i = j * 2;
+            t = GePrecomp::select(j, es[i]);
+            r = &h + &t;
+            h = r.to_full();
+        }
+
+        h
     }
 }
 
-impl Add<&GeCached> for &GeP3 {
+impl Add<&GeCached> for &Ge {
     type Output = GeP1P1;
 
     #[allow(clippy::suspicious_arithmetic_impl)]
@@ -269,15 +400,7 @@ impl Add<&GeCached> for &GeP3 {
     }
 }
 
-impl Add<GePrecomp> for GeP3 {
-    type Output = GeP1P1;
-
-    fn add(self, rhs: GePrecomp) -> GeP1P1 {
-        &self + &rhs
-    }
-}
-
-impl Add<&GePrecomp> for &GeP3 {
+impl Add<&GePrecomp> for &Ge {
     type Output = GeP1P1;
 
     #[allow(clippy::suspicious_arithmetic_impl)]
@@ -302,7 +425,7 @@ impl Add<&GePrecomp> for &GeP3 {
     }
 }
 
-impl Sub<GeCached> for GeP3 {
+impl Sub<GeCached> for Ge {
     type Output = GeP1P1;
 
     fn sub(self, rhs: GeCached) -> GeP1P1 {
@@ -310,7 +433,7 @@ impl Sub<GeCached> for GeP3 {
     }
 }
 
-impl Sub<&GeCached> for &GeP3 {
+impl Sub<&GeCached> for &Ge {
     type Output = GeP1P1;
 
     #[allow(clippy::suspicious_arithmetic_impl)]
@@ -336,7 +459,7 @@ impl Sub<&GeCached> for &GeP3 {
     }
 }
 
-impl Sub<GePrecomp> for GeP3 {
+impl Sub<GePrecomp> for Ge {
     type Output = GeP1P1;
 
     fn sub(self, rhs: GePrecomp) -> GeP1P1 {
@@ -344,7 +467,7 @@ impl Sub<GePrecomp> for GeP3 {
     }
 }
 
-impl Sub<&GePrecomp> for &GeP3 {
+impl Sub<&GePrecomp> for &Ge {
     type Output = GeP1P1;
 
     #[allow(clippy::suspicious_arithmetic_impl)]
@@ -370,13 +493,11 @@ impl Sub<&GePrecomp> for &GeP3 {
 }
 
 impl GePrecomp {
-    fn zero() -> GePrecomp {
-        GePrecomp {
-            y_plus_x: Fe::ONE,
-            y_minus_x: Fe::ONE,
-            xy2d: Fe::ZERO,
-        }
-    }
+    pub const ZERO: Self = Self {
+        y_plus_x: Fe::ONE,
+        y_minus_x: Fe::ONE,
+        xy2d: Fe::ZERO,
+    };
 
     pub(crate) fn maybe_set(&mut self, other: &GePrecomp, do_swap: Choice) {
         self.y_plus_x.maybe_set(&other.y_plus_x, do_swap);
@@ -389,7 +510,7 @@ impl GePrecomp {
 
         let bnegative = (b as u8) >> 7;
         let babs: u8 = (b - (((-(bnegative as i8)) & b) << 1)) as u8;
-        let mut t = GePrecomp::zero();
+        let mut t = GePrecomp::ZERO;
         t.maybe_set(&precomp::GE_BASE[pos][0], babs.ct_eq(1));
         t.maybe_set(&precomp::GE_BASE[pos][1], babs.ct_eq(2));
         t.maybe_set(&precomp::GE_BASE[pos][2], babs.ct_eq(3));
diff --git a/src/curve25519/mod.rs b/src/curve25519/mod.rs
index ba42fbf..afaf3f2 100644
--- a/src/curve25519/mod.rs
+++ b/src/curve25519/mod.rs
@@ -31,70 +31,11 @@ mod ge;
 pub mod scalar;
 
 pub use fe::Fe;
-pub use ge::{GeCached, GeP1P1, GeP2, GeP3, GePrecomp};
+pub use ge::{Ge, GeCached, GeP1P1, GePartial, GePrecomp};
 pub use scalar::Scalar;
 
 use crate::constant_time::CtZero;
 
-/*
-h = a * B
-where a = a[0]+256*a[1]+...+256^31 a[31]
-B is the Ed25519 base point (x,4/5) with x positive.
-
-Preconditions:
-  a[31] <= 127
-*/
-#[doc(hidden)]
-pub fn ge_scalarmult_base(a: &[u8; 32]) -> GeP3 {
-    let mut es: [i8; 64] = [0; 64];
-    let mut r: GeP1P1;
-    let mut s: GeP2;
-    let mut t: GePrecomp;
-
-    for i in 0..32 {
-        es[2 * i + 0] = ((a[i] >> 0) & 0b1111) as i8;
-        es[2 * i + 1] = ((a[i] >> 4) & 0b1111) as i8;
-    }
-    /* each es[i] is between 0 and 0xf */
-    /* es[63] is between 0 and 7 */
-
-    let mut carry: i8 = 0;
-    for esi in es[0..63].iter_mut() {
-        *esi += carry;
-        carry = *esi + 8;
-        carry >>= 4;
-        *esi -= carry << 4;
-    }
-    es[63] += carry;
-    /* each es[i] is between -8 and 8 */
-
-    let mut h = GeP3::ZERO;
-    for j in 0..32 {
-        let i = j * 2 + 1;
-        t = GePrecomp::select(j, es[i]);
-        r = h + t;
-        h = r.to_p3();
-    }
-
-    r = h.dbl();
-    s = r.to_p2();
-    r = s.dbl();
-    s = r.to_p2();
-    r = s.dbl();
-    s = r.to_p2();
-    r = s.dbl();
-    h = r.to_p3();
-
-    for j in 0..32 {
-        let i = j * 2;
-        t = GePrecomp::select(j, es[i]);
-        r = h + t;
-        h = r.to_p3();
-    }
-
-    h
-}
-
 /// Computes a shared secret from the curve25519 private key (n) and public
 /// key (p)
 pub fn curve25519(n: &[u8; 32], p: &[u8; 32]) -> [u8; 32] {
diff --git a/src/curve25519/scalar/mod.rs b/src/curve25519/scalar/mod.rs
index 8681514..a07de0e 100644
--- a/src/curve25519/scalar/mod.rs
+++ b/src/curve25519/scalar/mod.rs
@@ -136,14 +136,12 @@ mod tests {
         ];
 
         for (i, iv) in ivs.iter().enumerate() {
-            let mut out = [0u8; 32];
-            muladd(
-                &mut out,
+            let out = muladd(
                 &Scalar::from_bytes(&iv.a),
-                &iv.b,
+                &Scalar::from_bytes(&iv.b),
                 &Scalar::from_bytes(&iv.c),
             );
-            assert_eq!(iv.r, out, "IV test {} failed", i);
+            assert_eq!(iv.r, out.to_bytes(), "IV test {} failed", i);
         }
     }
 }
diff --git a/src/curve25519/scalar/scalar32.rs b/src/curve25519/scalar/scalar32.rs
index ec4253d..212232b 100644
--- a/src/curve25519/scalar/scalar32.rs
+++ b/src/curve25519/scalar/scalar32.rs
@@ -59,6 +59,20 @@ impl Scalar {
         r
     }
 
+    /// Get the scalar in a form of 64 nibbles
+    ///
+    /// nibble is a group of 4-bits
+    pub(crate) fn nibbles(&self) -> [i8; 64] {
+        let mut es: [i8; 64] = [0; 64];
+        let a = self.0;
+
+        for i in 0..32 {
+            es[2 * i + 0] = ((a[i] >> 0) & 0b1111) as i8;
+            es[2 * i + 1] = ((a[i] >> 4) & 0b1111) as i8;
+        }
+        es
+    }
+
     /// Create a new scalar from 64 bytes (512 bits) reducing
     /// the scalar to an element of the field
     ///
@@ -323,7 +337,7 @@ Output:
     where l = 2^252 + 27742317777372353535851937790883648493.
 */
 #[rustfmt::skip]
-pub(crate) fn muladd(s: &mut [u8; 32], Scalar(a): &Scalar, b: &[u8; 32], Scalar(c): &Scalar) {
+pub(crate) fn muladd(Scalar(a): &Scalar, Scalar(b): &Scalar, Scalar(c): &Scalar) -> Scalar {
     let a0 = 2097151 & load_3i(&a[0..3]);
     let a1 = 2097151 & (load_4i(&a[2..6]) >> 5);
     let a2 = 2097151 & (load_3i(&a[5..8]) >> 2);
@@ -622,6 +636,7 @@ pub(crate) fn muladd(s: &mut [u8; 32], Scalar(a): &Scalar, b: &[u8; 32], Scalar(
     carry9 = s9 >> 21; s10 += carry9; s9 -= carry9 << 21;
     carry10 = s10 >> 21; s11 += carry10; s10 -= carry10 << 21;
 
+    let mut s = [0u8; 32];
     s[0] = (s0 >> 0) as u8;
     s[1] = (s0 >> 8) as u8;
     s[2] = ((s0 >> 16) | (s1 << 5)) as u8;
@@ -654,4 +669,5 @@ pub(crate) fn muladd(s: &mut [u8; 32], Scalar(a): &Scalar, b: &[u8; 32], Scalar(
     s[29] = (s11 >> 1) as u8;
     s[30] = (s11 >> 9) as u8;
     s[31] = (s11 >> 17) as u8;
+    Scalar(s)
 }
diff --git a/src/curve25519/scalar/scalar64.rs b/src/curve25519/scalar/scalar64.rs
index d943420..3639743 100644
--- a/src/curve25519/scalar/scalar64.rs
+++ b/src/curve25519/scalar/scalar64.rs
@@ -161,6 +161,41 @@ impl Scalar {
         }
         r
     }
+
+    /// Get the scalar in a form of 64 nibbles
+    ///
+    /// nibble is a group of 4-bits
+    pub(crate) fn nibbles(&self) -> [i8; 64] {
+        let mut es: [i8; 64] = [0; 64];
+
+        // contract limbs
+        let mut c = [0; 4];
+        c[0] = self.0[1] << 56 | self.0[0];
+        c[1] = self.0[2] << 48 | self.0[1] >> 8;
+        c[2] = self.0[3] << 40 | self.0[2] >> 16;
+        c[3] = self.0[4] << 32 | self.0[3] >> 24;
+
+        // write 16 nibbles for each saturated limbs, for 64 nibbles
+        for b in 0..4 {
+            es[16 * b + 0] = ((c[b] >> 0) & 0b1111) as i8;
+            es[16 * b + 1] = ((c[b] >> 4) & 0b1111) as i8;
+            es[16 * b + 2] = ((c[b] >> 8) & 0b1111) as i8;
+            es[16 * b + 3] = ((c[b] >> 12) & 0b1111) as i8;
+            es[16 * b + 4] = ((c[b] >> 16) & 0b1111) as i8;
+            es[16 * b + 5] = ((c[b] >> 20) & 0b1111) as i8;
+            es[16 * b + 6] = ((c[b] >> 24) & 0b1111) as i8;
+            es[16 * b + 7] = ((c[b] >> 28) & 0b1111) as i8;
+            es[16 * b + 8] = ((c[b] >> 32) & 0b1111) as i8;
+            es[16 * b + 9] = ((c[b] >> 36) & 0b1111) as i8;
+            es[16 * b + 10] = ((c[b] >> 40) & 0b1111) as i8;
+            es[16 * b + 11] = ((c[b] >> 44) & 0b1111) as i8;
+            es[16 * b + 12] = ((c[b] >> 48) & 0b1111) as i8;
+            es[16 * b + 13] = ((c[b] >> 52) & 0b1111) as i8;
+            es[16 * b + 14] = ((c[b] >> 56) & 0b1111) as i8;
+            es[16 * b + 15] = ((c[b] >> 60) & 0b1111) as i8;
+        }
+        es
+    }
 }
 
 #[rustfmt::skip]
@@ -321,13 +356,10 @@ Output:
     where l = 2^252 + 27742317777372353535851937790883648493.
 */
 #[rustfmt::skip]
-pub(crate) fn muladd(s: &mut [u8; 32], a: &Scalar, b: &[u8; 32], c: &Scalar) {
-    let b = Scalar::from_bytes(b);
-
+pub(crate) fn muladd(a: &Scalar, b: &Scalar, c: &Scalar) -> Scalar {
     let m = mul(&a, &b);
     let r = add(&m, &c);
-
-    *s = r.to_bytes();
+    r
 }
 
 #[cfg(test)]
diff --git a/src/ed25519.rs b/src/ed25519.rs
index 00af886..2716dcc 100644
--- a/src/ed25519.rs
+++ b/src/ed25519.rs
@@ -22,9 +22,8 @@
 //!
 
 use crate::constant_time::CtEqual;
-use crate::curve25519::{curve25519, ge_scalarmult_base, scalar, Fe, GeP2, GeP3, Scalar};
-use crate::digest::Digest;
-use crate::sha2::Sha512;
+use crate::curve25519::{curve25519, scalar, Fe, Ge, GePartial, Scalar};
+use crate::hashing::sha2::Sha512;
 use core::convert::TryFrom;
 
 #[deprecated(since = "0.4.0", note = "use `PRIVATE_KEY_LENGTH`")]
@@ -50,10 +49,7 @@ fn clamp_scalar(scalar: &mut [u8]) {
 ///
 /// SCALAR(32bytes) | RANDOM(32bytes) = CLAMP(SHA512(private_key))
 fn extended_secret(private_key: &[u8; PRIVATE_KEY_LENGTH]) -> [u8; EXTENDED_KEY_LENGTH] {
-    let mut hash_output = [0; 64];
-    let mut hasher = Sha512::new();
-    hasher.input(private_key);
-    hasher.result(&mut hash_output);
+    let mut hash_output = Sha512::new().update(private_key).finalize();
     clamp_scalar(&mut hash_output);
     hash_output
 }
@@ -69,13 +65,18 @@ pub fn keypair_public(keypair: &[u8; KEYPAIR_LENGTH]) -> &[u8; PUBLIC_KEY_LENGTH
 }
 
 /// Extract the scalar part (first 32 bytes) from the extended key
-fn extended_scalar(extended_secret: &[u8; EXTENDED_KEY_LENGTH]) -> &[u8; 32] {
+fn extended_scalar(extended_secret: &[u8; EXTENDED_KEY_LENGTH]) -> Scalar {
+    Scalar::from_bytes(<&[u8; 32]>::try_from(&extended_secret[0..32]).unwrap())
+}
+
+/// Extract the scalar part (first 32 bytes) from the extended key
+fn extended_scalar_bytes(extended_secret: &[u8; EXTENDED_KEY_LENGTH]) -> &[u8; 32] {
     <&[u8; 32]>::try_from(&extended_secret[0..32]).unwrap()
 }
 
 /// generate the public key associated with an extended secret key
 pub fn extended_to_public(extended_secret: &[u8; EXTENDED_KEY_LENGTH]) -> [u8; PUBLIC_KEY_LENGTH] {
-    let a = ge_scalarmult_base(extended_scalar(extended_secret));
+    let a = Ge::scalarmult_base(&extended_scalar(extended_secret));
     a.to_bytes()
 }
 
@@ -100,11 +101,10 @@ pub fn keypair(
 /// Generate the nonce which is a scalar out of the extended_secret random part and the message itself
 /// using SHA512 and scalar_reduction
 fn signature_nonce(extended_secret: &[u8; EXTENDED_KEY_LENGTH], message: &[u8]) -> Scalar {
-    let mut hash_output: [u8; 64] = [0; 64];
-    let mut hasher = Sha512::new();
-    hasher.input(&extended_secret[32..64]);
-    hasher.input(message);
-    hasher.result(&mut hash_output);
+    let hash_output = Sha512::new()
+        .update(&extended_secret[32..64])
+        .update(message)
+        .finalize();
     Scalar::reduce_from_wide_bytes(&hash_output)
 }
 
@@ -116,25 +116,17 @@ pub fn signature(message: &[u8], keypair: &[u8; KEYPAIR_LENGTH]) -> [u8; SIGNATU
 
     let nonce = signature_nonce(&az, message);
 
-    let r: GeP3 = ge_scalarmult_base(&nonce.to_bytes());
+    let r = Ge::scalarmult_base(&nonce);
 
     let mut signature = [0; SIGNATURE_LENGTH];
     signature[0..32].copy_from_slice(&r.to_bytes());
     signature[32..64].copy_from_slice(public_key);
 
     {
-        let mut hasher = Sha512::new();
-        hasher.input(&signature);
-        hasher.input(message);
-        let mut hram: [u8; 64] = [0; 64];
-        hasher.result(&mut hram);
+        let hram = Sha512::new().update(&signature).update(message).finalize();
         let hram = Scalar::reduce_from_wide_bytes(&hram);
-        scalar::muladd(
-            <&mut [u8; 32]>::try_from(&mut signature[32..64]).unwrap(),
-            &hram,
-            extended_scalar(&az),
-            &nonce,
-        );
+        let r = scalar::muladd(&hram, &extended_scalar(&az), &nonce);
+        signature[32..64].copy_from_slice(&r.to_bytes())
     }
 
     signature
@@ -148,25 +140,17 @@ pub fn signature_extended(
     let public_key = extended_to_public(extended_secret);
     let nonce = signature_nonce(extended_secret, message);
 
-    let r: GeP3 = ge_scalarmult_base(&nonce.to_bytes());
+    let r = Ge::scalarmult_base(&nonce);
 
     let mut signature = [0; SIGNATURE_LENGTH];
     signature[0..32].copy_from_slice(&r.to_bytes());
     signature[32..64].copy_from_slice(&public_key);
 
     {
-        let mut hasher = Sha512::new();
-        hasher.input(&signature);
-        hasher.input(message);
-        let mut hram: [u8; 64] = [0; 64];
-        hasher.result(&mut hram);
-        let hram = Scalar::reduce_from_wide_bytes(&mut hram);
-        scalar::muladd(
-            <&mut [u8; 32]>::try_from(&mut signature[32..64]).unwrap(),
-            &hram,
-            extended_scalar(extended_secret),
-            &nonce,
-        );
+        let hram = Sha512::new().update(&signature).update(message).finalize();
+        let hram = Scalar::reduce_from_wide_bytes(&hram);
+        let r = scalar::muladd(&hram, &extended_scalar(extended_secret), &nonce);
+        signature[32..64].copy_from_slice(&r.to_bytes())
     }
 
     signature
@@ -181,7 +165,7 @@ pub fn verify(
     let signature_left = <&[u8; 32]>::try_from(&signature[0..32]).unwrap();
     let signature_right = <&[u8; 32]>::try_from(&signature[32..64]).unwrap();
 
-    let a = match GeP3::from_bytes_negate_vartime(public_key) {
+    let a = match Ge::from_bytes(public_key) {
         Some(g) => g,
         None => {
             return false;
@@ -193,6 +177,7 @@ pub fn verify(
         Some(s) => s,
     };
 
+    // reject all-0 public keys
     let mut d = 0;
     for pk_byte in public_key.iter() {
         d |= *pk_byte;
@@ -201,15 +186,14 @@ pub fn verify(
         return false;
     }
 
-    let mut hasher = Sha512::new();
-    hasher.input(signature_left);
-    hasher.input(public_key);
-    hasher.input(message);
-    let mut hash: [u8; 64] = [0; 64];
-    hasher.result(&mut hash);
-    let a_scalar = Scalar::reduce_from_wide_bytes(&mut hash);
+    let hash = Sha512::new()
+        .update(signature_left)
+        .update(public_key)
+        .update(message)
+        .finalize();
+    let a_scalar = Scalar::reduce_from_wide_bytes(&hash);
 
-    let r = GeP2::double_scalarmult_vartime(&a_scalar, a, &signature_scalar);
+    let r = GePartial::double_scalarmult_vartime(&a_scalar, a, &signature_scalar);
     let rcheck = r.to_bytes();
 
     CtEqual::ct_eq(&rcheck, signature_left).into()
@@ -224,7 +208,7 @@ pub fn exchange(public_key: &[u8; 32], private_key: &[u8; PRIVATE_KEY_LENGTH]) -
     // Produce private key from seed component (bytes 0 to 32)
     // of the Ed25519 extended private key (64 bytes).
     let extended_secret = extended_secret(private_key);
-    let shared_mont_x = curve25519(extended_scalar(&extended_secret), &mont_x.to_bytes());
+    let shared_mont_x = curve25519(extended_scalar_bytes(&extended_secret), &mont_x.to_bytes());
 
     shared_mont_x
 }