diff --git a/src/cryptographic_primitives/secret_sharing/feldman_vss.rs b/src/cryptographic_primitives/secret_sharing/feldman_vss.rs
index 9b6c29fc..014a13b0 100644
--- a/src/cryptographic_primitives/secret_sharing/feldman_vss.rs
+++ b/src/cryptographic_primitives/secret_sharing/feldman_vss.rs
@@ -242,37 +242,16 @@ impl<E: Curve> VerifiableSS<E> {
     //compute \lambda_{index,S}, a lagrangian coefficient that change the (t,n) scheme to (|S|,|S|)
     // used in http://stevengoldfeder.com/papers/GG18.pdf
     pub fn map_share_to_new_params(
-        params: &ShamirSecretSharing,
+        _params: &ShamirSecretSharing,
         index: u16,
         s: &[u16],
     ) -> Scalar<E> {
-        let s_len = s.len();
-        //     assert!(s_len > self.reconstruct_limit());
-        // add one to indices to get points
-        let points: Vec<Scalar<E>> = (0..params.share_count)
-            .map(|i| Scalar::from(i + 1))
-            .collect();
-
-        let xi = &points[usize::from(index)];
-        let num = Scalar::from(1);
-        let denum = Scalar::from(1);
-        let num = (0..s_len).fold(num, |acc, i| {
-            if s[i] != index {
-                acc * &points[usize::from(s[i])]
-            } else {
-                acc
-            }
-        });
-        let denum = (0..s_len).fold(denum, |acc, i| {
-            if s[i] != index {
-                let xj_sub_xi = &points[usize::from(s[i])] - xi;
-                acc * xj_sub_xi
-            } else {
-                acc
-            }
-        });
-        let denum = denum.invert().unwrap();
-        num * denum
+        let j = (0u16..)
+            .zip(s)
+            .find_map(|(j, s_j)| if *s_j == index { Some(j) } else { None })
+            .expect("`s` doesn't include `index`");
+        let xs = s.iter().map(|x| Scalar::from(*x + 1)).collect::<Vec<_>>();
+        Polynomial::lagrange_basis(&Scalar::zero(), j, &xs)
     }
 }
 
diff --git a/src/cryptographic_primitives/secret_sharing/polynomial.rs b/src/cryptographic_primitives/secret_sharing/polynomial.rs
index 96da82ed..9128f11c 100644
--- a/src/cryptographic_primitives/secret_sharing/polynomial.rs
+++ b/src/cryptographic_primitives/secret_sharing/polynomial.rs
@@ -230,6 +230,55 @@ impl<E: Curve> Polynomial<E> {
     pub fn coefficients(&self) -> &[Scalar<E>] {
         &self.coefficients
     }
+
+    /// Evaluates lagrange basis polynomial `l_j(x, xs)`
+    ///
+    /// Lagrange basis polynomials are mainly used for Lagrange interpolation, ie. calculating `L(x)`
+    /// where polynomial `L` is defined as set of `t+1` distinct points `(x_i, y_i)` (`t = deg(f)`).
+    /// [Example] section shows how Lagrange interpolation can be implemented using this function.
+    ///
+    /// [Example]: Self::lagrange_basis#example
+    ///
+    /// ## Panics
+    /// This function will panic if elements in `xs` are not pairwise distinct, or `j >= xs.len()`
+    ///
+    /// ## Example
+    /// If you have polynomial `f` defined as `t+1` points `(x_0, y_0), ..., (x_t, y_t)` (and polynomial
+    /// degree is `t`), then you can, for instance, calculate `f(15)`:
+    ///
+    /// ```rust
+    /// use curv::cryptographic_primitives::secret_sharing::Polynomial;
+    /// # use curv::elliptic::curves::*;
+    ///
+    /// # let t = 3;
+    /// # let f = Polynomial::<Secp256r1>::sample_exact(t);
+    /// # let (x_0, x_1, x_2, x_3) = (Scalar::from(1), Scalar::from(2), Scalar::from(3), Scalar::from(4));
+    /// # let (y_0, y_1, y_2, y_3) = (f.evaluate(&x_0), f.evaluate(&x_1), f.evaluate(&x_2), f.evaluate(&x_3));
+    /// let xs = &[x_0, x_1, x_2, x_3];
+    /// let ys = &[y_0, y_1, y_2, y_3];
+    ///
+    /// let f_15: Scalar<_> = (0..).zip(ys)
+    ///     .map(|(j, y_j)| y_j * Polynomial::lagrange_basis(&Scalar::from(15), j, xs))
+    ///     .sum();
+    /// assert_eq!(f_15, f.evaluate(&Scalar::from(15)));
+    /// ```
+    pub fn lagrange_basis(x: &Scalar<E>, j: u16, xs: &[Scalar<E>]) -> Scalar<E> {
+        let x_j = &xs[usize::from(j)];
+        let num: Scalar<E> = (0u16..)
+            .zip(xs)
+            .filter(|(m, _)| *m != j)
+            .map(|(_, x_m)| x - x_m)
+            .product();
+        let denum: Scalar<E> = (0u16..)
+            .zip(xs)
+            .filter(|(m, _)| *m != j)
+            .map(|(_, x_m)| x_j - x_m)
+            .product();
+        let denum = denum
+            .invert()
+            .expect("elements in xs are not pairwise distinct");
+        num * denum
+    }
 }
 
 /// Multiplies polynomial `f(x)` at scalar `s`, returning resulting polynomial `g(x) = s * f(x)`
diff --git a/src/elliptic/curves/p256.rs b/src/elliptic/curves/p256.rs
index 9ad20a9b..b93974a4 100644
--- a/src/elliptic/curves/p256.rs
+++ b/src/elliptic/curves/p256.rs
@@ -180,12 +180,15 @@ impl ECScalar for Secp256r1Scalar {
     }
 
     fn add_assign(&mut self, other: &Self) {
+        self.purpose = "add_assign";
         *self.fe += &*other.fe
     }
     fn mul_assign(&mut self, other: &Self) {
+        self.purpose = "mul_assign";
         *self.fe *= &*other.fe
     }
     fn sub_assign(&mut self, other: &Self) {
+        self.purpose = "sub_assign";
         *self.fe -= &*other.fe
     }