diff --git a/std/algebra/emulated/sw_emulated/point.go b/std/algebra/emulated/sw_emulated/point.go
index f09b34c553..329d64c5b4 100644
--- a/std/algebra/emulated/sw_emulated/point.go
+++ b/std/algebra/emulated/sw_emulated/point.go
@@ -604,14 +604,12 @@ func (c *Curve[B, S]) scalarMulGeneric(p *AffinePoint[B], s *emulated.Element[S]
 	if err != nil {
 		panic(fmt.Sprintf("parse opts: %v", err))
 	}
-	addFn := c.Add
 	var selector frontend.Variable
 	if cfg.UseSafe {
 		// if p=(0,0) we assign a dummy (0,1) to p and continue
 		selector = c.api.And(c.baseApi.IsZero(&p.X), c.baseApi.IsZero(&p.Y))
 		one := c.baseApi.One()
 		p = c.Select(selector, &AffinePoint[B]{X: *one, Y: *one}, p)
-		addFn = c.AddUnified
 	}
 
 	var st S
@@ -638,9 +636,10 @@ func (c *Curve[B, S]) scalarMulGeneric(p *AffinePoint[B], s *emulated.Element[S]
 	R0 = c.Select(sBits[n-1], Rb, R0)
 
 	// i = 0
-	// When cfg.UseSafe is set, we use AddUnified instead of Add. This means
-	// when s=0 then Acc=(0,0) because AddUnified(Q, -Q) = (0,0).
-	R0 = c.Select(sBits[0], R0, addFn(R0, c.Neg(p)))
+	// we use AddUnified instead of Add. This is because:
+	// 		- when s=0 then R0=P and AddUnified(P, -P) = (0,0). We return (0,0).
+	// 		- when s=1 then R0=P AddUnified(Q, -Q) is well defined. We return R0=P.
+	R0 = c.Select(sBits[0], R0, c.AddUnified(R0, c.Neg(p)))
 
 	if cfg.UseSafe {
 		// if p=(0,0), return (0,0)
@@ -665,26 +664,36 @@ func (c *Curve[B, S]) jointScalarMul(p1, p2 *AffinePoint[B], s1, s2 *emulated.El
 	}
 }
 
-// jointScalarMulGeneric computes s1 * p1 + s2 * p2 and returns it. It doesn't modify the inputs.
+// jointScalarMulGeneric computes [s1]p1 + [s2]p2. It doesn't modify the inputs.
 //
-// ⚠️  The scalar s must be nonzero and the point Q different from (0,0).
+// ⚠️  p1, p2 must not be (0,0) and s1, s2 must not be 0, unless [algopts.WithUseSafe] option is set.
 func (c *Curve[B, S]) jointScalarMulGeneric(p1, p2 *AffinePoint[B], s1, s2 *emulated.Element[S], opts ...algopts.AlgebraOption) *AffinePoint[B] {
-	s1r := c.scalarApi.Reduce(s1)
-	s1Bits := c.scalarApi.ToBits(s1r)
-	s2r := c.scalarApi.Reduce(s2)
-	s2Bits := c.scalarApi.ToBits(s2r)
-
-	res := c.scalarBitsMul(p1, s1Bits, opts...)
-	tmp := c.scalarBitsMul(p2, s2Bits, opts...)
-
-	res = c.Add(res, tmp)
-	return res
+	res1 := c.scalarMulGeneric(p1, s1, opts...)
+	res2 := c.scalarMulGeneric(p2, s2, opts...)
+	return c.Add(res1, res2)
 }
 
-// jointScalarMulGLV computes P = [s]Q + [t]R using Shamir's trick with an efficient endomorphism and returns it. It doesn't modify P, Q nor s.
+// jointScalarMulGLV computes [s1]p1 + [s2]p2 using an endomorphism. It doesn't modify P, Q nor s.
 //
-// ⚠️  The scalar s must be nonzero and the point Q different from (0,0).
-func (c *Curve[B, S]) jointScalarMulGLV(Q, R *AffinePoint[B], s, t *emulated.Element[S], opts ...algopts.AlgebraOption) *AffinePoint[B] {
+// ⚠️  The scalars s1, s2 must be nonzero and the point p1, p2 different from (0,0), unless [algopts.WithUseSafe] option is set.
+func (c *Curve[B, S]) jointScalarMulGLV(p1, p2 *AffinePoint[B], s1, s2 *emulated.Element[S], opts ...algopts.AlgebraOption) *AffinePoint[B] {
+	cfg, err := algopts.NewConfig(opts...)
+	if err != nil {
+		panic(fmt.Sprintf("parse opts: %v", err))
+	}
+	if cfg.UseSafe {
+		// TODO @yelhousni: optimize
+		res1 := c.scalarMulGLV(p1, s1, opts...)
+		res2 := c.scalarMulGLV(p2, s2, opts...)
+		return c.AddUnified(res1, res2)
+	} else {
+		return c.jointScalarMulGLVUnsafe(p1, p2, s1, s2)
+	}
+}
+
+// jointScalarMulGLVUnsafe computes [s]Q + [t]R using Shamir's trick with an efficient endomorphism and returns it. It doesn't modify P, Q nor s.
+// ⚠️  The scalar s must be nonzero and the point Q different from (0,0), unless [algopts.WithUseSafe] option is set.
+func (c *Curve[B, S]) jointScalarMulGLVUnsafe(Q, R *AffinePoint[B], s, t *emulated.Element[S]) *AffinePoint[B] {
 	var st S
 	frModulus := c.scalarApi.Modulus()
 	sd, err := c.scalarApi.NewHint(decomposeScalarG1, 5, s, c.eigenvalue, frModulus)
@@ -806,15 +815,23 @@ func (c *Curve[B, S]) jointScalarMulGLV(Q, R *AffinePoint[B], s, t *emulated.Ele
 	Acc = c.Select(t2bits[0], Acc, tablePhiR[0])
 
 	return Acc
+
 }
 
-// scalarBitsMul computes s * p and returns it where sBits is the bit decomposition of s. It doesn't modify p nor sBits.
-// ⚠️  Point and scalar must be nonzero.
-func (c *Curve[B, S]) scalarBitsMul(p *AffinePoint[B], sBits []frontend.Variable, opts ...algopts.AlgebraOption) *AffinePoint[B] {
+// scalarBitsMulGeneric computes s * p and returns it where sBits is the bit decomposition of s. It doesn't modify p nor sBits.
+// ⚠️  p must not be (0,0) and sBits not [0,...,0], unless [algopts.WithUseSafe] option is set.
+func (c *Curve[B, S]) scalarBitsMulGeneric(p *AffinePoint[B], sBits []frontend.Variable, opts ...algopts.AlgebraOption) *AffinePoint[B] {
 	cfg, err := algopts.NewConfig(opts...)
 	if err != nil {
 		panic(fmt.Sprintf("parse opts: %v", err))
 	}
+	var selector frontend.Variable
+	if cfg.UseSafe {
+		// if p=(0,0) we assign a dummy (0,1) to p and continue
+		selector = c.api.And(c.baseApi.IsZero(&p.X), c.baseApi.IsZero(&p.Y))
+		one := c.baseApi.One()
+		p = c.Select(selector, &AffinePoint[B]{X: *one, Y: *one}, p)
+	}
 
 	var st S
 	n := st.Modulus().BitLen()
@@ -838,8 +855,17 @@ func (c *Curve[B, S]) scalarBitsMul(p *AffinePoint[B], sBits []frontend.Variable
 	R0 = c.Select(sBits[n-1], Rb, R0)
 
 	// i = 0
+	// we use AddUnified instead of Add. This is because:
+	// 		- when s=0 then R0=P and AddUnified(P, -P) = (0,0). We return (0,0).
+	// 		- when s=1 then R0=P AddUnified(Q, -Q) is well defined. We return R0=P.
 	R0 = c.Select(sBits[0], R0, c.AddUnified(R0, c.Neg(p)))
 
+	if cfg.UseSafe {
+		// if p=(0,0), return (0,0)
+		zero := c.baseApi.Zero()
+		R0 = c.Select(selector, &AffinePoint[B]{X: *zero, Y: *zero}, R0)
+	}
+
 	return R0
 }
 
@@ -884,7 +910,13 @@ func (c *Curve[B, S]) ScalarMulBase(s *emulated.Element[S], opts ...algopts.Alge
 	}
 
 	// i = 0
-	tmp := c.AddUnified(res, c.Neg(g))
+	// When cfg.UseSafe is set, we use AddUnified instead of Add. This means
+	// when s=0 then Acc=(0,0) because AddUnified(Q, -Q) = (0,0).
+	addFn := c.Add
+	if cfg.UseSafe {
+		addFn = c.AddUnified
+	}
+	tmp := addFn(res, c.Neg(g))
 	res = c.Select(sBits[0], res, tmp)
 
 	return res
@@ -911,8 +943,8 @@ func (c *Curve[B, S]) ScalarMulBase(s *emulated.Element[S], opts ...algopts.Alge
 //
 // This saves the Select logic related to (0,0) and the use of AddUnified to
 // handle the 0-scalar edge case.
-func (c *Curve[B, S]) JointScalarMulBase(p *AffinePoint[B], s2, s1 *emulated.Element[S]) *AffinePoint[B] {
-	return c.jointScalarMul(c.Generator(), p, s1, s2)
+func (c *Curve[B, S]) JointScalarMulBase(p *AffinePoint[B], s2, s1 *emulated.Element[S], opts ...algopts.AlgebraOption) *AffinePoint[B] {
+	return c.jointScalarMul(c.Generator(), p, s1, s2, opts...)
 }
 
 // MultiScalarMul computes the multi scalar multiplication of the points P and
@@ -932,6 +964,10 @@ func (c *Curve[B, S]) MultiScalarMul(p []*AffinePoint[B], s []*emulated.Element[
 	if err != nil {
 		return nil, fmt.Errorf("new config: %w", err)
 	}
+	addFn := c.Add
+	if cfg.UseSafe {
+		addFn = c.AddUnified
+	}
 	if !cfg.FoldMulti {
 		// the scalars are unique
 		if len(p) != len(s) {
@@ -946,7 +982,7 @@ func (c *Curve[B, S]) MultiScalarMul(p []*AffinePoint[B], s []*emulated.Element[
 		}
 		for i := 1; i < n-1; i += 2 {
 			q := c.jointScalarMul(p[i-1], p[i], s[i-1], s[i], opts...)
-			res = c.Add(res, q)
+			res = addFn(res, q)
 		}
 		return res, nil
 	} else {
@@ -957,12 +993,12 @@ func (c *Curve[B, S]) MultiScalarMul(p []*AffinePoint[B], s []*emulated.Element[
 		gamma := s[0]
 		gamma = c.scalarApi.Reduce(gamma)
 		gammaBits := c.scalarApi.ToBits(gamma)
-		res := c.scalarBitsMul(p[len(p)-1], gammaBits, opts...)
+		res := c.scalarBitsMulGeneric(p[len(p)-1], gammaBits, opts...)
 		for i := len(p) - 2; i > 0; i-- {
-			res = c.Add(p[i], res)
-			res = c.scalarBitsMul(res, gammaBits, opts...)
+			res = addFn(p[i], res)
+			res = c.scalarBitsMulGeneric(res, gammaBits, opts...)
 		}
-		res = c.Add(p[0], res)
+		res = addFn(p[0], res)
 		return res, nil
 	}
 }
diff --git a/std/algebra/emulated/sw_emulated/point_test.go b/std/algebra/emulated/sw_emulated/point_test.go
index d6dba62ddc..859e63dcb8 100644
--- a/std/algebra/emulated/sw_emulated/point_test.go
+++ b/std/algebra/emulated/sw_emulated/point_test.go
@@ -923,6 +923,137 @@ func TestJointScalarMulBase(t *testing.T) {
 	assert.NoError(err)
 }
 
+type MultiScalarMulEdgeCasesTest[T, S emulated.FieldParams] struct {
+	Points  []AffinePoint[T]
+	Scalars []emulated.Element[S]
+	Res     AffinePoint[T]
+}
+
+func (c *MultiScalarMulEdgeCasesTest[T, S]) Define(api frontend.API) error {
+	cr, err := New[T, S](api, GetCurveParams[T]())
+	if err != nil {
+		return err
+	}
+	ps := make([]*AffinePoint[T], len(c.Points))
+	for i := range c.Points {
+		ps[i] = &c.Points[i]
+	}
+	ss := make([]*emulated.Element[S], len(c.Scalars))
+	for i := range c.Scalars {
+		ss[i] = &c.Scalars[i]
+	}
+	res, err := cr.MultiScalarMul(ps, ss, algopts.WithUseSafe())
+	if err != nil {
+		return err
+	}
+	cr.AssertIsEqual(res, &c.Res)
+	return nil
+}
+
+func TestMultiScalarMulEdgeCases(t *testing.T) {
+	assert := test.NewAssert(t)
+	nbLen := 5
+	P := make([]bw6761.G1Affine, nbLen)
+	S := make([]fr_bw6761.Element, nbLen)
+	for i := 0; i < nbLen; i++ {
+		S[i].SetRandom()
+		P[i].ScalarMultiplicationBase(S[i].BigInt(new(big.Int)))
+	}
+	var res bw6761.G1Affine
+	_, err := res.MultiExp(P, S, ecc.MultiExpConfig{})
+
+	assert.NoError(err)
+	cP := make([]AffinePoint[emulated.BW6761Fp], len(P))
+	cS := make([]emulated.Element[emparams.BW6761Fr], len(S))
+	var infinity bw6761.G1Affine
+
+	// s1 * (0,0) + s2 * (0,0) + s3 * (0,0) + s4 * (0,0)  + s5 * (0,0) == (0,0)
+	for i := range cP {
+		cP[i] = AffinePoint[emparams.BW6761Fp]{
+			X: emulated.ValueOf[emparams.BW6761Fp](infinity.X),
+			Y: emulated.ValueOf[emparams.BW6761Fp](infinity.Y),
+		}
+	}
+	for i := range cS {
+		cS[i] = emulated.ValueOf[emparams.BW6761Fr](S[i])
+	}
+	assignment1 := MultiScalarMulEdgeCasesTest[emparams.BW6761Fp, emparams.BW6761Fr]{
+		Points:  cP,
+		Scalars: cS,
+		Res: AffinePoint[emparams.BW6761Fp]{
+			X: emulated.ValueOf[emparams.BW6761Fp](infinity.X),
+			Y: emulated.ValueOf[emparams.BW6761Fp](infinity.Y),
+		},
+	}
+	err = test.IsSolved(&MultiScalarMulEdgeCasesTest[emparams.BW6761Fp, emparams.BW6761Fr]{
+		Points:  make([]AffinePoint[emparams.BW6761Fp], nbLen),
+		Scalars: make([]emulated.Element[emparams.BW6761Fr], nbLen),
+	}, &assignment1, ecc.BN254.ScalarField())
+	assert.NoError(err)
+
+	// 0 * P1 + 0 * P2 + 0 * P3 + 0 * P4 + 0 * P5 == (0,0)
+	for i := range cP {
+		cP[i] = AffinePoint[emparams.BW6761Fp]{
+			X: emulated.ValueOf[emparams.BW6761Fp](P[i].X),
+			Y: emulated.ValueOf[emparams.BW6761Fp](P[i].Y),
+		}
+	}
+	for i := range cS {
+		cS[i] = emulated.ValueOf[emparams.BW6761Fr](0)
+	}
+	assignment2 := MultiScalarMulEdgeCasesTest[emparams.BW6761Fp, emparams.BW6761Fr]{
+		Points:  cP,
+		Scalars: cS,
+		Res: AffinePoint[emparams.BW6761Fp]{
+			X: emulated.ValueOf[emparams.BW6761Fp](infinity.X),
+			Y: emulated.ValueOf[emparams.BW6761Fp](infinity.Y),
+		},
+	}
+	err = test.IsSolved(&MultiScalarMulEdgeCasesTest[emparams.BW6761Fp, emparams.BW6761Fr]{
+		Points:  make([]AffinePoint[emparams.BW6761Fp], nbLen),
+		Scalars: make([]emulated.Element[emparams.BW6761Fr], nbLen),
+	}, &assignment2, ecc.BN254.ScalarField())
+	assert.NoError(err)
+
+	// s1 * (0,0) + s2 * P2 + s3 * (0,0) + s4 * P4 + 0 * P5 == s2 * P + s4 * P4
+	var res3 bw6761.G1Affine
+	res3.ScalarMultiplication(&P[1], S[1].BigInt(new(big.Int)))
+	res.ScalarMultiplication(&P[3], S[3].BigInt(new(big.Int)))
+	res3.Add(&res3, &res)
+	for i := range cP {
+		cP[i] = AffinePoint[emparams.BW6761Fp]{
+			X: emulated.ValueOf[emparams.BW6761Fp](P[i].X),
+			Y: emulated.ValueOf[emparams.BW6761Fp](P[i].Y),
+		}
+	}
+	cP[0] = AffinePoint[emparams.BW6761Fp]{
+		X: emulated.ValueOf[emparams.BW6761Fp](infinity.X),
+		Y: emulated.ValueOf[emparams.BW6761Fp](infinity.Y),
+	}
+	cP[2] = AffinePoint[emparams.BW6761Fp]{
+		X: emulated.ValueOf[emparams.BW6761Fp](infinity.X),
+		Y: emulated.ValueOf[emparams.BW6761Fp](infinity.Y),
+	}
+	for i := range cS {
+		cS[i] = emulated.ValueOf[emparams.BW6761Fr](S[i])
+	}
+	cS[4] = emulated.ValueOf[emparams.BW6761Fr](0)
+	assignment3 := MultiScalarMulEdgeCasesTest[emparams.BW6761Fp, emparams.BW6761Fr]{
+		Points:  cP,
+		Scalars: cS,
+		Res: AffinePoint[emparams.BW6761Fp]{
+			X: emulated.ValueOf[emparams.BW6761Fp](res3.X),
+			Y: emulated.ValueOf[emparams.BW6761Fp](res3.Y),
+		},
+	}
+	err = test.IsSolved(&MultiScalarMulEdgeCasesTest[emparams.BW6761Fp, emparams.BW6761Fr]{
+		Points:  make([]AffinePoint[emparams.BW6761Fp], nbLen),
+		Scalars: make([]emulated.Element[emparams.BW6761Fr], nbLen),
+	}, &assignment3, ecc.BN254.ScalarField())
+	assert.NoError(err)
+
+}
+
 type MultiScalarMulTest[T, S emulated.FieldParams] struct {
 	Points  []AffinePoint[T]
 	Scalars []emulated.Element[S]
@@ -989,6 +1120,176 @@ func TestMultiScalarMul(t *testing.T) {
 	assert.NoError(err)
 }
 
+type MultiScalarMulFoldedEdgeCasesTest[T, S emulated.FieldParams] struct {
+	Points  []AffinePoint[T]
+	Scalars []emulated.Element[S]
+	Res     AffinePoint[T]
+}
+
+func (c *MultiScalarMulFoldedEdgeCasesTest[T, S]) Define(api frontend.API) error {
+	cr, err := New[T, S](api, GetCurveParams[T]())
+	if err != nil {
+		return err
+	}
+	ps := make([]*AffinePoint[T], len(c.Points))
+	for i := range c.Points {
+		ps[i] = &c.Points[i]
+	}
+	ss := make([]*emulated.Element[S], len(c.Scalars))
+	for i := range c.Scalars {
+		ss[i] = &c.Scalars[i]
+	}
+	res, err := cr.MultiScalarMul(ps, ss, algopts.WithFoldingScalarMul(), algopts.WithUseSafe())
+	if err != nil {
+		return err
+	}
+	cr.AssertIsEqual(res, &c.Res)
+	return nil
+}
+
+func TestMultiScalarFoldedEdgeCasesMul(t *testing.T) {
+	assert := test.NewAssert(t)
+	nbLen := 5
+	P := make([]bw6761.G1Affine, nbLen)
+	S := make([]fr_bw6761.Element, nbLen)
+	S[0].SetOne()
+	S[1].SetRandom()
+	S[2].Square(&S[1])
+	S[3].Mul(&S[1], &S[2])
+	S[4].Mul(&S[1], &S[3])
+	for i := 0; i < nbLen; i++ {
+		P[i].ScalarMultiplicationBase(S[i].BigInt(new(big.Int)))
+	}
+	var res, infinity bw6761.G1Affine
+	_, err := res.MultiExp(P, S, ecc.MultiExpConfig{})
+
+	assert.NoError(err)
+	cP := make([]AffinePoint[emulated.BW6761Fp], len(P))
+	cS := make([]emulated.Element[emparams.BW6761Fr], len(S))
+
+	// s^0 * (0,0) + s^1 * (0,0) + s^2 * (0,0) + s^3 * (0,0)  + s^4 * (0,0) == (0,0)
+	for i := range cP {
+		cP[i] = AffinePoint[emparams.BW6761Fp]{
+			X: emulated.ValueOf[emparams.BW6761Fp](infinity.X),
+			Y: emulated.ValueOf[emparams.BW6761Fp](infinity.Y),
+		}
+	}
+	// s0 = s
+	S[0].Set(&S[1])
+	for i := range cS {
+		cS[i] = emulated.ValueOf[emparams.BW6761Fr](S[i])
+	}
+	assignment1 := MultiScalarMulFoldedEdgeCasesTest[emparams.BW6761Fp, emparams.BW6761Fr]{
+		Points:  cP,
+		Scalars: cS,
+		Res: AffinePoint[emparams.BW6761Fp]{
+			X: emulated.ValueOf[emparams.BW6761Fp](infinity.X),
+			Y: emulated.ValueOf[emparams.BW6761Fp](infinity.Y),
+		},
+	}
+	err = test.IsSolved(&MultiScalarMulFoldedEdgeCasesTest[emparams.BW6761Fp, emparams.BW6761Fr]{
+		Points:  make([]AffinePoint[emparams.BW6761Fp], nbLen),
+		Scalars: make([]emulated.Element[emparams.BW6761Fr], nbLen),
+	}, &assignment1, ecc.BN254.ScalarField())
+	assert.NoError(err)
+
+	// 0^0 * P1 + 0 * P2 + 0 * P3 + 0 * P4 + 0 * P5 == P1
+	for i := range cP {
+		cP[i] = AffinePoint[emparams.BW6761Fp]{
+			X: emulated.ValueOf[emparams.BW6761Fp](P[i].X),
+			Y: emulated.ValueOf[emparams.BW6761Fp](P[i].Y),
+		}
+	}
+	for i := range cS {
+		cS[i] = emulated.ValueOf[emparams.BW6761Fr](0)
+	}
+	assignment2 := MultiScalarMulFoldedEdgeCasesTest[emparams.BW6761Fp, emparams.BW6761Fr]{
+		Points:  cP,
+		Scalars: cS,
+		Res: AffinePoint[emparams.BW6761Fp]{
+			X: emulated.ValueOf[emparams.BW6761Fp](P[0].X),
+			Y: emulated.ValueOf[emparams.BW6761Fp](P[0].Y),
+		},
+	}
+	err = test.IsSolved(&MultiScalarMulFoldedEdgeCasesTest[emparams.BW6761Fp, emparams.BW6761Fr]{
+		Points:  make([]AffinePoint[emparams.BW6761Fp], nbLen),
+		Scalars: make([]emulated.Element[emparams.BW6761Fr], nbLen),
+	}, &assignment2, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+type MultiScalarMulFoldedTest[T, S emulated.FieldParams] struct {
+	Points  []AffinePoint[T]
+	Scalars []emulated.Element[S]
+	Res     AffinePoint[T]
+}
+
+func (c *MultiScalarMulFoldedTest[T, S]) Define(api frontend.API) error {
+	cr, err := New[T, S](api, GetCurveParams[T]())
+	if err != nil {
+		return err
+	}
+	ps := make([]*AffinePoint[T], len(c.Points))
+	for i := range c.Points {
+		ps[i] = &c.Points[i]
+	}
+	ss := make([]*emulated.Element[S], len(c.Scalars))
+	for i := range c.Scalars {
+		ss[i] = &c.Scalars[i]
+	}
+	res, err := cr.MultiScalarMul(ps, ss, algopts.WithFoldingScalarMul())
+	if err != nil {
+		return err
+	}
+	cr.AssertIsEqual(res, &c.Res)
+	return nil
+}
+
+func TestMultiScalarFoldedMul(t *testing.T) {
+	assert := test.NewAssert(t)
+	nbLen := 4
+	P := make([]bw6761.G1Affine, nbLen)
+	S := make([]fr_bw6761.Element, nbLen)
+	// [s^0]P0 + [s^1]P1 + [s^2]P2 + [s^3]P3 = P0 + [s]P1 + [s^2]P2 + [s^3]P3
+	S[0].SetOne()
+	S[1].SetRandom()
+	S[2].Square(&S[1])
+	S[3].Mul(&S[1], &S[2])
+	for i := 0; i < nbLen; i++ {
+		P[i].ScalarMultiplicationBase(S[i].BigInt(new(big.Int)))
+	}
+	var res bw6761.G1Affine
+	_, err := res.MultiExp(P, S, ecc.MultiExpConfig{})
+
+	assert.NoError(err)
+	cP := make([]AffinePoint[emulated.BW6761Fp], len(P))
+	for i := range cP {
+		cP[i] = AffinePoint[emparams.BW6761Fp]{
+			X: emulated.ValueOf[emparams.BW6761Fp](P[i].X),
+			Y: emulated.ValueOf[emparams.BW6761Fp](P[i].Y),
+		}
+	}
+	cS := make([]emulated.Element[emparams.BW6761Fr], len(S))
+	// s0 = s
+	S[0].Set(&S[1])
+	for i := range cS {
+		cS[i] = emulated.ValueOf[emparams.BW6761Fr](S[i])
+	}
+	assignment := MultiScalarMulFoldedTest[emparams.BW6761Fp, emparams.BW6761Fr]{
+		Points:  cP,
+		Scalars: cS,
+		Res: AffinePoint[emparams.BW6761Fp]{
+			X: emulated.ValueOf[emparams.BW6761Fp](res.X),
+			Y: emulated.ValueOf[emparams.BW6761Fp](res.Y),
+		},
+	}
+	err = test.IsSolved(&MultiScalarMulFoldedTest[emparams.BW6761Fp, emparams.BW6761Fr]{
+		Points:  make([]AffinePoint[emparams.BW6761Fp], nbLen),
+		Scalars: make([]emulated.Element[emparams.BW6761Fr], nbLen),
+	}, &assignment, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
 type ScalarMulTestBounded[T, S emulated.FieldParams] struct {
 	P, Q AffinePoint[T]
 	S    emulated.Element[S]
@@ -1038,8 +1339,6 @@ func TestScalarMulBounded(t *testing.T) {
 	assert.NoError(err)
 }
 
-//
-
 type JointScalarMulTest[T, S emulated.FieldParams] struct {
 	P1, P2, Q AffinePoint[T]
 	S1, S2    emulated.Element[S]
@@ -1092,6 +1391,159 @@ func TestJointScalarMul6(t *testing.T) {
 	assert.NoError(err)
 }
 
+type JointScalarMulEdgeCasesTest[T, S emulated.FieldParams] struct {
+	P1, P2, Q AffinePoint[T]
+	S1, S2    emulated.Element[S]
+}
+
+func (c *JointScalarMulEdgeCasesTest[T, S]) Define(api frontend.API) error {
+	cr, err := New[T, S](api, GetCurveParams[T]())
+	if err != nil {
+		return err
+	}
+	res := cr.jointScalarMul(&c.P1, &c.P2, &c.S1, &c.S2, algopts.WithUseSafe())
+	cr.AssertIsEqual(res, &c.Q)
+	return nil
+}
+
+func TestJointScalarMulEdgeCases6(t *testing.T) {
+	assert := test.NewAssert(t)
+	var r1, r2 fr_bw6761.Element
+	_, _ = r1.SetRandom()
+	_, _ = r2.SetRandom()
+	s1 := new(big.Int)
+	s2 := new(big.Int)
+	r1.BigInt(s1)
+	r2.BigInt(s2)
+	var res, res1, res2, gen2, infinity bw6761.G1Affine
+	_, _, gen1, _ := bw6761.Generators()
+	gen2.Double(&gen1)
+	res1.ScalarMultiplication(&gen1, s1)
+	res2.ScalarMultiplication(&gen2, s2)
+	res.Add(&res1, &res2)
+
+	circuit := JointScalarMulEdgeCasesTest[emulated.BW6761Fp, emulated.BW6761Fr]{}
+	// s1*(0,0) + s2*(0,0) == (0,0)
+	witness1 := JointScalarMulTest[emulated.BW6761Fp, emulated.BW6761Fr]{
+		S1: emulated.ValueOf[emulated.BW6761Fr](s1),
+		S2: emulated.ValueOf[emulated.BW6761Fr](s2),
+		P1: AffinePoint[emulated.BW6761Fp]{
+			X: emulated.ValueOf[emulated.BW6761Fp](infinity.X),
+			Y: emulated.ValueOf[emulated.BW6761Fp](infinity.Y),
+		},
+		P2: AffinePoint[emulated.BW6761Fp]{
+			X: emulated.ValueOf[emulated.BW6761Fp](infinity.X),
+			Y: emulated.ValueOf[emulated.BW6761Fp](infinity.Y),
+		},
+		Q: AffinePoint[emulated.BW6761Fp]{
+			X: emulated.ValueOf[emulated.BW6761Fp](infinity.X),
+			Y: emulated.ValueOf[emulated.BW6761Fp](infinity.Y),
+		},
+	}
+	err := test.IsSolved(&circuit, &witness1, testCurve.ScalarField())
+	assert.NoError(err)
+
+	// s1*P + s2*(0,0) == s1*P
+	witness2 := JointScalarMulTest[emulated.BW6761Fp, emulated.BW6761Fr]{
+		S1: emulated.ValueOf[emulated.BW6761Fr](s1),
+		S2: emulated.ValueOf[emulated.BW6761Fr](s2),
+		P1: AffinePoint[emulated.BW6761Fp]{
+			X: emulated.ValueOf[emulated.BW6761Fp](gen1.X),
+			Y: emulated.ValueOf[emulated.BW6761Fp](gen1.Y),
+		},
+		P2: AffinePoint[emulated.BW6761Fp]{
+			X: emulated.ValueOf[emulated.BW6761Fp](infinity.X),
+			Y: emulated.ValueOf[emulated.BW6761Fp](infinity.Y),
+		},
+		Q: AffinePoint[emulated.BW6761Fp]{
+			X: emulated.ValueOf[emulated.BW6761Fp](res1.X),
+			Y: emulated.ValueOf[emulated.BW6761Fp](res1.Y),
+		},
+	}
+	err = test.IsSolved(&circuit, &witness2, testCurve.ScalarField())
+	assert.NoError(err)
+
+	// s1*(0,0) + s2*Q == s2*Q
+	witness3 := JointScalarMulTest[emulated.BW6761Fp, emulated.BW6761Fr]{
+		S1: emulated.ValueOf[emulated.BW6761Fr](s1),
+		S2: emulated.ValueOf[emulated.BW6761Fr](s2),
+		P1: AffinePoint[emulated.BW6761Fp]{
+			X: emulated.ValueOf[emulated.BW6761Fp](infinity.X),
+			Y: emulated.ValueOf[emulated.BW6761Fp](infinity.Y),
+		},
+		P2: AffinePoint[emulated.BW6761Fp]{
+			X: emulated.ValueOf[emulated.BW6761Fp](gen2.X),
+			Y: emulated.ValueOf[emulated.BW6761Fp](gen2.Y),
+		},
+		Q: AffinePoint[emulated.BW6761Fp]{
+			X: emulated.ValueOf[emulated.BW6761Fp](res2.X),
+			Y: emulated.ValueOf[emulated.BW6761Fp](res2.Y),
+		},
+	}
+	err = test.IsSolved(&circuit, &witness3, testCurve.ScalarField())
+	assert.NoError(err)
+
+	// 0*P + 0*Q == (0,0)
+	witness4 := JointScalarMulTest[emulated.BW6761Fp, emulated.BW6761Fr]{
+		S1: emulated.ValueOf[emulated.BW6761Fr](0),
+		S2: emulated.ValueOf[emulated.BW6761Fr](0),
+		P1: AffinePoint[emulated.BW6761Fp]{
+			X: emulated.ValueOf[emulated.BW6761Fp](gen1.X),
+			Y: emulated.ValueOf[emulated.BW6761Fp](gen1.Y),
+		},
+		P2: AffinePoint[emulated.BW6761Fp]{
+			X: emulated.ValueOf[emulated.BW6761Fp](gen2.X),
+			Y: emulated.ValueOf[emulated.BW6761Fp](gen2.Y),
+		},
+		Q: AffinePoint[emulated.BW6761Fp]{
+			X: emulated.ValueOf[emulated.BW6761Fp](infinity.X),
+			Y: emulated.ValueOf[emulated.BW6761Fp](infinity.Y),
+		},
+	}
+	err = test.IsSolved(&circuit, &witness4, testCurve.ScalarField())
+	assert.NoError(err)
+
+	// 0*P + s2*Q == s2*Q
+	witness5 := JointScalarMulTest[emulated.BW6761Fp, emulated.BW6761Fr]{
+		S1: emulated.ValueOf[emulated.BW6761Fr](0),
+		S2: emulated.ValueOf[emulated.BW6761Fr](s2),
+		P1: AffinePoint[emulated.BW6761Fp]{
+			X: emulated.ValueOf[emulated.BW6761Fp](gen1.X),
+			Y: emulated.ValueOf[emulated.BW6761Fp](gen1.Y),
+		},
+		P2: AffinePoint[emulated.BW6761Fp]{
+			X: emulated.ValueOf[emulated.BW6761Fp](gen2.X),
+			Y: emulated.ValueOf[emulated.BW6761Fp](gen2.Y),
+		},
+		Q: AffinePoint[emulated.BW6761Fp]{
+			X: emulated.ValueOf[emulated.BW6761Fp](res2.X),
+			Y: emulated.ValueOf[emulated.BW6761Fp](res2.Y),
+		},
+	}
+	err = test.IsSolved(&circuit, &witness5, testCurve.ScalarField())
+	assert.NoError(err)
+
+	// s1*P + 0*Q == s1*P
+	witness6 := JointScalarMulTest[emulated.BW6761Fp, emulated.BW6761Fr]{
+		S1: emulated.ValueOf[emulated.BW6761Fr](s1),
+		S2: emulated.ValueOf[emulated.BW6761Fr](0),
+		P1: AffinePoint[emulated.BW6761Fp]{
+			X: emulated.ValueOf[emulated.BW6761Fp](gen1.X),
+			Y: emulated.ValueOf[emulated.BW6761Fp](gen1.Y),
+		},
+		P2: AffinePoint[emulated.BW6761Fp]{
+			X: emulated.ValueOf[emulated.BW6761Fp](gen2.X),
+			Y: emulated.ValueOf[emulated.BW6761Fp](gen2.Y),
+		},
+		Q: AffinePoint[emulated.BW6761Fp]{
+			X: emulated.ValueOf[emulated.BW6761Fp](res1.X),
+			Y: emulated.ValueOf[emulated.BW6761Fp](res1.Y),
+		},
+	}
+	err = test.IsSolved(&circuit, &witness6, testCurve.ScalarField())
+	assert.NoError(err)
+}
+
 type MuxCircuitTest[T, S emulated.FieldParams] struct {
 	Selector frontend.Variable
 	Inputs   [8]AffinePoint[T]
diff --git a/std/algebra/native/sw_bls12377/g1.go b/std/algebra/native/sw_bls12377/g1.go
index 6040a297ea..d484959cc0 100644
--- a/std/algebra/native/sw_bls12377/g1.go
+++ b/std/algebra/native/sw_bls12377/g1.go
@@ -28,24 +28,11 @@ import (
 	"github.com/consensys/gnark/std/algebra/algopts"
 )
 
-// G1Jac point in Jacobian coords
-type G1Jac struct {
-	X, Y, Z frontend.Variable
-}
-
 // G1Affine point in affine coords
 type G1Affine struct {
 	X, Y frontend.Variable
 }
 
-// Neg outputs -p
-func (p *G1Jac) Neg(api frontend.API, p1 G1Jac) *G1Jac {
-	p.X = p1.X
-	p.Y = api.Sub(0, p1.Y)
-	p.Z = p1.Z
-	return p
-}
-
 // Neg outputs -p
 func (p *G1Affine) Neg(api frontend.API, p1 G1Affine) *G1Affine {
 	p.X = p1.X
@@ -113,91 +100,6 @@ func (p *G1Affine) AddUnified(api frontend.API, q G1Affine) *G1Affine {
 	return p
 }
 
-// AddAssign adds 2 point in Jacobian coordinates
-// p=p, a=p1
-func (p *G1Jac) AddAssign(api frontend.API, p1 G1Jac) *G1Jac {
-
-	// get some Element from our pool
-	var Z1Z1, Z2Z2, U1, U2, S1, S2, H, I, J, r, V frontend.Variable
-
-	Z1Z1 = api.Mul(p1.Z, p1.Z)
-
-	Z2Z2 = api.Mul(p.Z, p.Z)
-
-	U1 = api.Mul(p1.X, Z2Z2)
-
-	U2 = api.Mul(p.X, Z1Z1)
-
-	S1 = api.Mul(p1.Y, p.Z)
-	S1 = api.Mul(S1, Z2Z2)
-
-	S2 = api.Mul(p.Y, p1.Z)
-	S2 = api.Mul(S2, Z1Z1)
-
-	H = api.Sub(U2, U1)
-
-	I = api.Add(H, H)
-	I = api.Mul(I, I)
-
-	J = api.Mul(H, I)
-
-	r = api.Sub(S2, S1)
-	r = api.Add(r, r)
-
-	V = api.Mul(U1, I)
-
-	p.X = api.Mul(r, r)
-	p.X = api.Sub(p.X, J)
-	p.X = api.Sub(p.X, V)
-	p.X = api.Sub(p.X, V)
-
-	p.Y = api.Sub(V, p.X)
-	p.Y = api.Mul(p.Y, r)
-
-	S1 = api.Mul(J, S1)
-	S1 = api.Add(S1, S1)
-
-	p.Y = api.Sub(p.Y, S1)
-
-	p.Z = api.Add(p.Z, p1.Z)
-	p.Z = api.Mul(p.Z, p.Z)
-	p.Z = api.Sub(p.Z, Z1Z1)
-	p.Z = api.Sub(p.Z, Z2Z2)
-	p.Z = api.Mul(p.Z, H)
-
-	return p
-}
-
-// DoubleAssign doubles the receiver point in jacobian coords and returns it
-func (p *G1Jac) DoubleAssign(api frontend.API) *G1Jac {
-	// get some Element from our pool
-	var XX, YY, YYYY, ZZ, S, M, T frontend.Variable
-
-	XX = api.Mul(p.X, p.X)
-	YY = api.Mul(p.Y, p.Y)
-	YYYY = api.Mul(YY, YY)
-	ZZ = api.Mul(p.Z, p.Z)
-	S = api.Add(p.X, YY)
-	S = api.Mul(S, S)
-	S = api.Sub(S, XX)
-	S = api.Sub(S, YYYY)
-	S = api.Add(S, S)
-	M = api.Mul(XX, 3) // M = 3*XX+a*ZZ², here a=0 (we suppose sw has j invariant 0)
-	p.Z = api.Add(p.Z, p.Y)
-	p.Z = api.Mul(p.Z, p.Z)
-	p.Z = api.Sub(p.Z, YY)
-	p.Z = api.Sub(p.Z, ZZ)
-	p.X = api.Mul(M, M)
-	T = api.Add(S, S)
-	p.X = api.Sub(p.X, T)
-	p.Y = api.Sub(S, p.X)
-	p.Y = api.Mul(p.Y, M)
-	YYYY = api.Mul(YYYY, 8)
-	p.Y = api.Sub(p.Y, YYYY)
-
-	return p
-}
-
 // Select sets p1 if b=1, p2 if b=0, and returns it. b must be boolean constrained
 func (p *G1Affine) Select(api frontend.API, b frontend.Variable, p1, p2 G1Affine) *G1Affine {
 
@@ -223,14 +125,6 @@ func (p *G1Affine) Lookup2(api frontend.API, b1, b2 frontend.Variable, p1, p2, p
 
 }
 
-// FromJac sets p to p1 in affine and returns it
-func (p *G1Affine) FromJac(api frontend.API, p1 G1Jac) *G1Affine {
-	s := api.Mul(p1.Z, p1.Z)
-	p.X = api.DivUnchecked(p1.X, s)
-	p.Y = api.DivUnchecked(p1.Y, api.Mul(s, p1.Z))
-	return p
-}
-
 // Double double a point in affine coords
 func (p *G1Affine) Double(api frontend.API, p1 G1Affine) *G1Affine {
 
@@ -258,11 +152,11 @@ func (p *G1Affine) Double(api frontend.API, p1 G1Affine) *G1Affine {
 // The method chooses an implementation based on scalar s. If it is constant,
 // then the compiled circuit depends on s. If it is variable type, then
 // the circuit is independent of the inputs.
-func (P *G1Affine) ScalarMul(api frontend.API, Q G1Affine, s interface{}) *G1Affine {
+func (P *G1Affine) ScalarMul(api frontend.API, Q G1Affine, s interface{}, opts ...algopts.AlgebraOption) *G1Affine {
 	if n, ok := api.Compiler().ConstantValue(s); ok {
-		return P.constScalarMul(api, Q, n)
+		return P.constScalarMul(api, Q, n, opts...)
 	} else {
-		return P.varScalarMul(api, Q, s)
+		return P.varScalarMul(api, Q, s, opts...)
 	}
 }
 
@@ -292,7 +186,11 @@ func init() {
 }
 
 // varScalarMul sets P = [s] Q and returns P.
-func (P *G1Affine) varScalarMul(api frontend.API, Q G1Affine, s frontend.Variable) *G1Affine {
+func (P *G1Affine) varScalarMul(api frontend.API, Q G1Affine, s frontend.Variable, opts ...algopts.AlgebraOption) *G1Affine {
+	cfg, err := algopts.NewConfig(opts...)
+	if err != nil {
+		panic(err)
+	}
 	// This method computes [s] Q. We use several methods to reduce the number
 	// of added constraints - first, instead of classical double-and-add, we use
 	// the optimized version from https://github.com/zcash/zcash/issues/3924
@@ -304,7 +202,12 @@ func (P *G1Affine) varScalarMul(api frontend.API, Q G1Affine, s frontend.Variabl
 	// from a precomputed table. However, precomputing the table adds 12
 	// additional constraints and thus table-version is more expensive than
 	// addition-version.
-
+	var selector frontend.Variable
+	if cfg.UseSafe {
+		// if Q=(0,0) we assign a dummy (1,1) to Q and continue
+		selector = api.And(api.IsZero(Q.X), api.IsZero(Q.Y))
+		Q.Select(api, selector, G1Affine{X: 1, Y: 1}, Q)
+	}
 	// The context we are working is based on the `outer` curve. However, the
 	// points and the operations on the points are performed on the `inner`
 	// curve of the outer curve. We require some parameters from the inner
@@ -384,10 +287,21 @@ func (P *G1Affine) varScalarMul(api frontend.API, Q G1Affine, s frontend.Variabl
 		Acc.DoubleAndAdd(api, &Acc, &B)
 	}
 
-	tableQ[0].AddAssign(api, Acc)
-	Acc.Select(api, s1bits[0], Acc, tableQ[0])
-	tablePhiQ[0].AddAssign(api, Acc)
-	Acc.Select(api, s2bits[0], Acc, tablePhiQ[0])
+	// i = 0
+	// When cfg.UseSafe is set, we use AddUnified instead of Add. This means
+	// when s=0 then Acc=(0,0) because AddUnified(Q, -Q) = (0,0).
+	if cfg.UseSafe {
+		tableQ[0].AddUnified(api, Acc)
+		Acc.Select(api, s1bits[0], Acc, tableQ[0])
+		tablePhiQ[0].AddUnified(api, Acc)
+		Acc.Select(api, s2bits[0], Acc, tablePhiQ[0])
+		Acc.Select(api, selector, G1Affine{X: 0, Y: 0}, Acc)
+	} else {
+		tableQ[0].AddAssign(api, Acc)
+		Acc.Select(api, s1bits[0], Acc, tableQ[0])
+		tablePhiQ[0].AddAssign(api, Acc)
+		Acc.Select(api, s2bits[0], Acc, tablePhiQ[0])
+	}
 
 	P.X = Acc.X
 	P.Y = Acc.Y
@@ -396,7 +310,16 @@ func (P *G1Affine) varScalarMul(api frontend.API, Q G1Affine, s frontend.Variabl
 }
 
 // constScalarMul sets P = [s] Q and returns P.
-func (P *G1Affine) constScalarMul(api frontend.API, Q G1Affine, s *big.Int) *G1Affine {
+func (P *G1Affine) constScalarMul(api frontend.API, Q G1Affine, s *big.Int, opts ...algopts.AlgebraOption) *G1Affine {
+	cfg, err := algopts.NewConfig(opts...)
+	if err != nil {
+		panic(err)
+	}
+	if s.BitLen() == 0 {
+		P.X = 0
+		P.Y = 0
+		return P
+	}
 	// see the comments in varScalarMul. However, two-bit lookup is cheaper if
 	// bits are constant and here it makes sense to use the table in the main
 	// loop.
@@ -422,49 +345,60 @@ func (P *G1Affine) constScalarMul(api frontend.API, Q G1Affine, s *big.Int) *G1A
 	negPhiQ.Neg(api, phiQ)
 	var table [4]G1Affine
 	table[0] = negQ
-	table[0].AddAssign(api, negPhiQ)
 	table[1] = Q
-	table[1].AddAssign(api, negPhiQ)
 	table[2] = negQ
-	table[2].AddAssign(api, phiQ)
 	table[3] = Q
-	table[3].AddAssign(api, phiQ)
+
+	if cfg.UseSafe {
+		table[0].AddUnified(api, negPhiQ)
+		table[1].AddUnified(api, negPhiQ)
+		table[2].AddUnified(api, phiQ)
+		table[3].AddUnified(api, phiQ)
+	} else {
+		table[0].AddAssign(api, negPhiQ)
+		table[1].AddAssign(api, negPhiQ)
+		table[2].AddAssign(api, phiQ)
+		table[3].AddAssign(api, phiQ)
+	}
 
 	Acc = table[3]
 	// if both high bits are set, then we would get to the incomplete part,
 	// handle it separately.
 	if k[0].Bit(nbits-1) == 1 && k[1].Bit(nbits-1) == 1 {
-		Acc.Double(api, Acc)
-		Acc.AddAssign(api, table[3])
+		if cfg.UseSafe {
+			Acc.AddUnified(api, Acc)
+			Acc.AddUnified(api, table[3])
+		} else {
+			Acc.Double(api, Acc)
+			Acc.AddAssign(api, table[3])
+		}
 		nbits = nbits - 1
 	}
 	for i := nbits - 1; i > 0; i-- {
-		Acc.DoubleAndAdd(api, &Acc, &table[k[0].Bit(i)+2*k[1].Bit(i)])
+		if cfg.UseSafe {
+			Acc.AddUnified(api, Acc)
+			Acc.AddUnified(api, table[k[0].Bit(i)+2*k[1].Bit(i)])
+		} else {
+			Acc.DoubleAndAdd(api, &Acc, &table[k[0].Bit(i)+2*k[1].Bit(i)])
+		}
 	}
 
-	negQ.AddAssign(api, Acc)
-	Acc.Select(api, k[0].Bit(0), Acc, negQ)
-	negPhiQ.AddAssign(api, Acc)
+	// i = 0
+	if cfg.UseSafe {
+		negQ.AddUnified(api, Acc)
+		Acc.Select(api, k[0].Bit(0), Acc, negQ)
+		negPhiQ.AddUnified(api, Acc)
+	} else {
+		negQ.AddAssign(api, Acc)
+		Acc.Select(api, k[0].Bit(0), Acc, negQ)
+		negPhiQ.AddAssign(api, Acc)
+	}
 	Acc.Select(api, k[1].Bit(0), Acc, negPhiQ)
 	P.X, P.Y = Acc.X, Acc.Y
 
 	return P
 }
 
-// Assign a value to self (witness assignment)
-func (p *G1Jac) Assign(p1 *bls12377.G1Jac) {
-	p.X = (fr.Element)(p1.X)
-	p.Y = (fr.Element)(p1.Y)
-	p.Z = (fr.Element)(p1.Z)
-}
-
-// AssertIsEqual constraint self to be equal to other into the given constraint system
-func (p *G1Jac) AssertIsEqual(api frontend.API, other G1Jac) {
-	api.AssertIsEqual(p.X, other.X)
-	api.AssertIsEqual(p.Y, other.Y)
-	api.AssertIsEqual(p.Z, other.Z)
-}
-
 // Assign a value to self (witness assignment)
 func (p *G1Affine) Assign(p1 *bls12377.G1Affine) {
 	p.X = (fr.Element)(p1.X)
@@ -510,17 +444,34 @@ func (p *G1Affine) DoubleAndAdd(api frontend.API, p1, p2 *G1Affine) *G1Affine {
 }
 
 // ScalarMulBase computes s * g1 and returns it, where g1 is the fixed generator. It doesn't modify s.
-func (P *G1Affine) ScalarMulBase(api frontend.API, s frontend.Variable) *G1Affine {
+func (P *G1Affine) ScalarMulBase(api frontend.API, s frontend.Variable, opts ...algopts.AlgebraOption) *G1Affine {
 	_, _, g1aff, _ := bls12377.Generators()
 	generator := G1Affine{
 		X: g1aff.X.BigInt(new(big.Int)),
 		Y: g1aff.Y.BigInt(new(big.Int)),
 	}
-	return P.ScalarMul(api, generator, s)
+	return P.ScalarMul(api, generator, s, opts...)
+}
+
+func (P *G1Affine) jointScalarMul(api frontend.API, Q, R G1Affine, s, t frontend.Variable, opts ...algopts.AlgebraOption) *G1Affine {
+	cfg, err := algopts.NewConfig(opts...)
+	if err != nil {
+		panic(err)
+	}
+	if cfg.UseSafe {
+		// TODO @yelhousni: optimize
+		var tmp G1Affine
+		P.ScalarMul(api, Q, s, opts...)
+		tmp.ScalarMul(api, R, t, opts...)
+		P.AddUnified(api, tmp)
+	} else {
+		P.jointScalarMulUnsafe(api, Q, R, s, t)
+	}
+	return P
 }
 
 // P = [s]Q + [t]R using Shamir's trick
-func (P *G1Affine) jointScalarMul(api frontend.API, Q, R G1Affine, s, t frontend.Variable) *G1Affine {
+func (P *G1Affine) jointScalarMulUnsafe(api frontend.API, Q, R G1Affine, s, t frontend.Variable) *G1Affine {
 	cc := getInnerCurveConfig(api.Compiler().Field())
 
 	sd, err := api.Compiler().NewHint(DecomposeScalarG1, 3, s)
@@ -605,8 +556,19 @@ func (P *G1Affine) jointScalarMul(api frontend.API, Q, R G1Affine, s, t frontend
 	return P
 }
 
-// scalarBitsMul...
+// scalarBitsMul computes s * p and returns it where sBits is the bit decomposition of s. It doesn't modify p nor sBits.
+// The method is similar to varScalarMul.
 func (P *G1Affine) scalarBitsMul(api frontend.API, Q G1Affine, s1bits, s2bits []frontend.Variable, opts ...algopts.AlgebraOption) *G1Affine {
+	cfg, err := algopts.NewConfig(opts...)
+	if err != nil {
+		panic(err)
+	}
+	var selector frontend.Variable
+	if cfg.UseSafe {
+		// if Q=(0,0) we assign a dummy (1,1) to Q and continue
+		selector = api.And(api.IsZero(Q.X), api.IsZero(Q.Y))
+		Q.Select(api, selector, G1Affine{X: 1, Y: 1}, Q)
+	}
 	cc := getInnerCurveConfig(api.Compiler().Field())
 	nbits := cc.lambda.BitLen() + 1
 	var Acc /*accumulator*/, B, B2 /*tmp vars*/ G1Affine
@@ -655,10 +617,21 @@ func (P *G1Affine) scalarBitsMul(api frontend.API, Q G1Affine, s1bits, s2bits []
 		Acc.DoubleAndAdd(api, &Acc, &B)
 	}
 
-	tableQ[0].AddAssign(api, Acc)
-	Acc.Select(api, s1bits[0], Acc, tableQ[0])
-	tablePhiQ[0].AddAssign(api, Acc)
-	Acc.Select(api, s2bits[0], Acc, tablePhiQ[0])
+	// i = 0
+	// When cfg.UseSafe is set, we use AddUnified instead of Add. This means
+	// when s=0 then Acc=(0,0) because AddUnified(Q, -Q) = (0,0).
+	if cfg.UseSafe {
+		tableQ[0].AddUnified(api, Acc)
+		Acc.Select(api, s1bits[0], Acc, tableQ[0])
+		tablePhiQ[0].AddUnified(api, Acc)
+		Acc.Select(api, s2bits[0], Acc, tablePhiQ[0])
+		Acc.Select(api, selector, G1Affine{X: 0, Y: 0}, Acc)
+	} else {
+		tableQ[0].AddAssign(api, Acc)
+		Acc.Select(api, s1bits[0], Acc, tableQ[0])
+		tablePhiQ[0].AddAssign(api, Acc)
+		Acc.Select(api, s2bits[0], Acc, tablePhiQ[0])
+	}
 
 	P.X = Acc.X
 	P.Y = Acc.Y
diff --git a/std/algebra/native/sw_bls12377/g1_test.go b/std/algebra/native/sw_bls12377/g1_test.go
index 01bfe73345..0d832351f0 100644
--- a/std/algebra/native/sw_bls12377/g1_test.go
+++ b/std/algebra/native/sw_bls12377/g1_test.go
@@ -24,7 +24,9 @@ import (
 	"github.com/consensys/gnark-crypto/ecc/bls12-377/fp"
 	"github.com/consensys/gnark-crypto/ecc/bls12-377/fr"
 	"github.com/consensys/gnark/frontend"
+	"github.com/consensys/gnark/std/algebra/algopts"
 	"github.com/consensys/gnark/std/math/emulated"
+	"github.com/consensys/gnark/std/math/emulated/emparams"
 	"github.com/consensys/gnark/test"
 
 	bls12377 "github.com/consensys/gnark-crypto/ecc/bls12-377"
@@ -116,43 +118,6 @@ func TestMarshalG1(t *testing.T) {
 	})
 }
 
-// -------------------------------------------------------------------------------------------------
-// Add jacobian
-
-type g1AddAssign struct {
-	A, B G1Jac
-	C    G1Jac `gnark:",public"`
-}
-
-func (circuit *g1AddAssign) Define(api frontend.API) error {
-	expected := circuit.A
-	expected.AddAssign(api, circuit.B)
-	expected.AssertIsEqual(api, circuit.C)
-	return nil
-}
-
-func TestAddAssignG1(t *testing.T) {
-
-	// sample 2 random points
-	a := randomPointG1()
-	b := randomPointG1()
-
-	// create the cs
-	var circuit, witness g1AddAssign
-
-	// assign the inputs
-	witness.A.Assign(&a)
-	witness.B.Assign(&b)
-
-	// compute the result
-	a.AddAssign(&b)
-	witness.C.Assign(&a)
-
-	assert := test.NewAssert(t)
-	assert.CheckCircuit(&circuit, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_761))
-
-}
-
 // -------------------------------------------------------------------------------------------------
 // Add affine
 
@@ -194,41 +159,6 @@ func TestAddAssignAffineG1(t *testing.T) {
 
 }
 
-// -------------------------------------------------------------------------------------------------
-// Double Jacobian
-
-type g1DoubleAssign struct {
-	A G1Jac
-	C G1Jac `gnark:",public"`
-}
-
-func (circuit *g1DoubleAssign) Define(api frontend.API) error {
-	expected := circuit.A
-	expected.DoubleAssign(api)
-	expected.AssertIsEqual(api, circuit.C)
-	return nil
-}
-
-func TestDoubleAssignG1(t *testing.T) {
-
-	// sample 2 random points
-	a := randomPointG1()
-
-	// create the cs
-	var circuit, witness g1DoubleAssign
-
-	// assign the inputs
-	witness.A.Assign(&a)
-
-	// compute the result
-	a.DoubleAssign()
-	witness.C.Assign(&a)
-
-	assert := test.NewAssert(t)
-	assert.CheckCircuit(&circuit, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_761))
-
-}
-
 // -------------------------------------------------------------------------------------------------
 // Double affine
 
@@ -306,60 +236,71 @@ func TestDoubleAndAddAffineG1(t *testing.T) {
 }
 
 // -------------------------------------------------------------------------------------------------
-// Neg
+// Scalar multiplication
 
-type g1Neg struct {
-	A G1Jac
-	C G1Jac `gnark:",public"`
+type g1constantScalarMul struct {
+	A G1Affine
+	C G1Affine `gnark:",public"`
+	R *big.Int
 }
 
-func (circuit *g1Neg) Define(api frontend.API) error {
-	expected := G1Jac{}
-	expected.Neg(api, circuit.A)
+func (circuit *g1constantScalarMul) Define(api frontend.API) error {
+	expected := G1Affine{}
+	expected.constScalarMul(api, circuit.A, circuit.R)
 	expected.AssertIsEqual(api, circuit.C)
 	return nil
 }
 
-func TestNegG1(t *testing.T) {
-
-	// sample 2 random points
-	a := randomPointG1()
+func TestConstantScalarMulG1(t *testing.T) {
+	// sample random point
+	_a := randomPointG1()
+	var a, c bls12377.G1Affine
+	a.FromJacobian(&_a)
 
+	// create the cs
+	var circuit, witness g1constantScalarMul
+	var r fr.Element
+	_, _ = r.SetRandom()
 	// assign the inputs
-	var witness g1Neg
 	witness.A.Assign(&a)
-	a.Neg(&a)
-	witness.C.Assign(&a)
+	// compute the result
+	br := new(big.Int)
+	r.BigInt(br)
+	// br is a circuit parameter
+	circuit.R = br
+	_a.ScalarMultiplication(&_a, br)
+	c.FromJacobian(&_a)
+	witness.C.Assign(&c)
 
 	assert := test.NewAssert(t)
-	assert.SolvingSucceeded(&g1Neg{}, &witness, test.WithCurves(ecc.BW6_761))
+	assert.CheckCircuit(&circuit, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_761))
 
 }
 
-// -------------------------------------------------------------------------------------------------
-// Scalar multiplication
-
-type g1constantScalarMul struct {
+type g1constantScalarMulEdgeCases struct {
 	A G1Affine
-	C G1Affine `gnark:",public"`
 	R *big.Int
 }
 
-func (circuit *g1constantScalarMul) Define(api frontend.API) error {
-	expected := G1Affine{}
-	expected.constScalarMul(api, circuit.A, circuit.R)
-	expected.AssertIsEqual(api, circuit.C)
+func (circuit *g1constantScalarMulEdgeCases) Define(api frontend.API) error {
+	expected1 := G1Affine{}
+	expected2 := G1Affine{}
+	infinity := G1Affine{X: 0, Y: 0}
+	expected1.constScalarMul(api, circuit.A, big.NewInt(0))
+	expected2.constScalarMul(api, infinity, circuit.R, algopts.WithUseSafe())
+	expected1.AssertIsEqual(api, infinity)
+	expected2.AssertIsEqual(api, infinity)
 	return nil
 }
 
-func TestConstantScalarMulG1(t *testing.T) {
+func TestConstantScalarMulG1EdgeCases(t *testing.T) {
 	// sample random point
 	_a := randomPointG1()
-	var a, c bls12377.G1Affine
+	var a bls12377.G1Affine
 	a.FromJacobian(&_a)
 
 	// create the cs
-	var circuit, witness g1constantScalarMul
+	var circuit, witness g1constantScalarMulEdgeCases
 	var r fr.Element
 	_, _ = r.SetRandom()
 	// assign the inputs
@@ -369,9 +310,6 @@ func TestConstantScalarMulG1(t *testing.T) {
 	r.BigInt(br)
 	// br is a circuit parameter
 	circuit.R = br
-	_a.ScalarMultiplication(&_a, br)
-	c.FromJacobian(&_a)
-	witness.C.Assign(&c)
 
 	assert := test.NewAssert(t)
 	assert.CheckCircuit(&circuit, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_761))
@@ -414,6 +352,40 @@ func TestVarScalarMulG1(t *testing.T) {
 	assert.CheckCircuit(&circuit, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_761))
 }
 
+type g1varScalarMulEdgeCases struct {
+	A G1Affine
+	R frontend.Variable
+}
+
+func (circuit *g1varScalarMulEdgeCases) Define(api frontend.API) error {
+	expected1 := G1Affine{}
+	expected2 := G1Affine{}
+	infinity := G1Affine{X: 0, Y: 0}
+	expected1.varScalarMul(api, circuit.A, 0, algopts.WithUseSafe())
+	expected2.varScalarMul(api, infinity, circuit.R, algopts.WithUseSafe())
+	expected1.AssertIsEqual(api, infinity)
+	expected2.AssertIsEqual(api, infinity)
+	return nil
+}
+
+func TestVarScalarMulG1EdgeCases(t *testing.T) {
+	// sample random point
+	_a := randomPointG1()
+	var a bls12377.G1Affine
+	a.FromJacobian(&_a)
+
+	// create the cs
+	var circuit, witness g1varScalarMulEdgeCases
+	var r fr.Element
+	_, _ = r.SetRandom()
+	witness.R = r.String()
+	// assign the inputs
+	witness.A.Assign(&a)
+
+	assert := test.NewAssert(t)
+	assert.CheckCircuit(&circuit, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_761))
+}
+
 type g1ScalarMul struct {
 	A    G1Affine
 	C    G1Affine `gnark:",public"`
@@ -485,6 +457,114 @@ func TestVarScalarMulBaseG1(t *testing.T) {
 	assert.CheckCircuit(&circuit, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_761))
 }
 
+type MultiScalarMulEdgeCasesTest struct {
+	Points  []G1Affine
+	Scalars []emulated.Element[ScalarField]
+	Res     G1Affine
+}
+
+func (c *MultiScalarMulEdgeCasesTest) Define(api frontend.API) error {
+	cr, err := NewCurve(api)
+	if err != nil {
+		return err
+	}
+	ps := make([]*G1Affine, len(c.Points))
+	for i := range c.Points {
+		ps[i] = &c.Points[i]
+	}
+	ss := make([]*emulated.Element[ScalarField], len(c.Scalars))
+	for i := range c.Scalars {
+		ss[i] = &c.Scalars[i]
+	}
+	res, err := cr.MultiScalarMul(ps, ss, algopts.WithUseSafe())
+	if err != nil {
+		return err
+	}
+	cr.AssertIsEqual(res, &c.Res)
+	return nil
+}
+
+func TestMultiScalarMulEdgeCases(t *testing.T) {
+	assert := test.NewAssert(t)
+	nbLen := 5
+	P := make([]bls12377.G1Affine, nbLen)
+	S := make([]fr.Element, nbLen)
+	for i := 0; i < nbLen; i++ {
+		S[i].SetRandom()
+		P[i].ScalarMultiplicationBase(S[i].BigInt(new(big.Int)))
+	}
+	var res, infinity bls12377.G1Affine
+	_, err := res.MultiExp(P, S, ecc.MultiExpConfig{})
+
+	assert.NoError(err)
+	cP := make([]G1Affine, len(P))
+	cS := make([]emulated.Element[ScalarField], len(S))
+
+	// s1 * (0,0) + s2 * (0,0) + s3 * (0,0) + s4 * (0,0)  + s5 * (0,0) == (0,0)
+	for i := range cP {
+		cP[i] = NewG1Affine(infinity)
+	}
+	for i := range cS {
+		cS[i] = NewScalar(S[i])
+	}
+	assignment1 := MultiScalarMulEdgeCasesTest{
+		Points:  cP,
+		Scalars: cS,
+		Res:     NewG1Affine(infinity),
+	}
+	err = test.IsSolved(&MultiScalarMulEdgeCasesTest{
+		Points:  make([]G1Affine, nbLen),
+		Scalars: make([]emulated.Element[ScalarField], nbLen),
+	}, &assignment1, ecc.BW6_761.ScalarField())
+	assert.NoError(err)
+
+	// 0 * P1 + 0 * P2 + 0 * P3 + 0 * P4 + 0 * P5 == (0,0)
+	for i := range cP {
+		cP[i] = NewG1Affine(P[i])
+	}
+	for i := range cS {
+		cS[i] = emulated.ValueOf[emparams.BLS12377Fr](0)
+	}
+	assignment2 := MultiScalarMulEdgeCasesTest{
+		Points:  cP,
+		Scalars: cS,
+		Res:     NewG1Affine(infinity),
+	}
+	err = test.IsSolved(&MultiScalarMulEdgeCasesTest{
+		Points:  make([]G1Affine, nbLen),
+		Scalars: make([]emulated.Element[ScalarField], nbLen),
+	}, &assignment2, ecc.BW6_761.ScalarField())
+	assert.NoError(err)
+
+	// s1 * (0,0) + s2 * P2 + s3 * (0,0) + s4 * P4 + 0 * P5 == s2 * P + s4 * P4
+	var res3 bls12377.G1Affine
+	res3.ScalarMultiplication(&P[1], S[1].BigInt(new(big.Int)))
+	res.ScalarMultiplication(&P[3], S[3].BigInt(new(big.Int)))
+	res3.Add(&res3, &res)
+	for i := range cP {
+		cP[i] = NewG1Affine(P[i])
+	}
+	cP[0].X = infinity.X
+	cP[0].Y = infinity.Y
+	cP[2].X = infinity.X
+	cP[2].Y = infinity.Y
+	for i := range cS {
+		cS[i] = NewScalar(S[i])
+	}
+	cS[4] = emulated.ValueOf[emparams.BLS12377Fr](0)
+
+	assignment3 := MultiScalarMulEdgeCasesTest{
+		Points:  cP,
+		Scalars: cS,
+		Res:     NewG1Affine(res3),
+	}
+	err = test.IsSolved(&MultiScalarMulEdgeCasesTest{
+		Points:  make([]G1Affine, nbLen),
+		Scalars: make([]emulated.Element[ScalarField], nbLen),
+	}, &assignment3, ecc.BW6_761.ScalarField())
+	assert.NoError(err)
+}
+
 type MultiScalarMulTest struct {
 	Points  []G1Affine
 	Scalars []emulated.Element[ScalarField]
@@ -545,6 +625,61 @@ func TestMultiScalarMul(t *testing.T) {
 	assert.NoError(err)
 }
 
+type g1JointScalarMulEdgeCases struct {
+	A, B G1Affine
+	C    G1Affine `gnark:",public"`
+	R, S frontend.Variable
+}
+
+func (circuit *g1JointScalarMulEdgeCases) Define(api frontend.API) error {
+	expected1 := G1Affine{}
+	expected2 := G1Affine{}
+	expected3 := G1Affine{}
+	expected4 := G1Affine{}
+	infinity := G1Affine{X: 0, Y: 0}
+	expected1.jointScalarMul(api, infinity, infinity, circuit.R, circuit.S, algopts.WithUseSafe())
+	expected2.jointScalarMul(api, circuit.A, circuit.B, big.NewInt(0), big.NewInt(0), algopts.WithUseSafe())
+	expected3.jointScalarMul(api, circuit.A, infinity, circuit.R, circuit.S, algopts.WithUseSafe())
+	expected4.jointScalarMul(api, circuit.A, circuit.B, circuit.R, big.NewInt(0), algopts.WithUseSafe())
+	_expected := G1Affine{}
+	_expected.ScalarMul(api, circuit.A, circuit.R, algopts.WithUseSafe())
+	expected1.AssertIsEqual(api, infinity)
+	expected2.AssertIsEqual(api, infinity)
+	expected3.AssertIsEqual(api, _expected)
+	expected4.AssertIsEqual(api, _expected)
+	return nil
+}
+
+func TestJointScalarMulG1EdgeCases(t *testing.T) {
+	// sample random point
+	_a := randomPointG1()
+	_b := randomPointG1()
+	var a, b, c bls12377.G1Affine
+	a.FromJacobian(&_a)
+	b.FromJacobian(&_b)
+
+	// create the cs
+	var circuit, witness g1JointScalarMulEdgeCases
+	var r, s fr.Element
+	_, _ = r.SetRandom()
+	_, _ = s.SetRandom()
+	witness.R = r.String()
+	witness.S = s.String()
+	// assign the inputs
+	witness.A.Assign(&a)
+	witness.B.Assign(&b)
+	// compute the result
+	var br, bs big.Int
+	_a.ScalarMultiplication(&_a, r.BigInt(&br))
+	_b.ScalarMultiplication(&_b, s.BigInt(&bs))
+	_a.AddAssign(&_b)
+	c.FromJacobian(&_a)
+	witness.C.Assign(&c)
+
+	assert := test.NewAssert(t)
+	assert.CheckCircuit(&circuit, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_761))
+}
+
 type g1JointScalarMul struct {
 	A, B G1Affine
 	C    G1Affine `gnark:",public"`
@@ -645,3 +780,156 @@ func randomPointG1() bls12377.G1Jac {
 
 	return p1
 }
+
+type MultiScalarMulFoldedEdgeCasesTest struct {
+	Points  []G1Affine
+	Scalars []emulated.Element[ScalarField]
+	Res     G1Affine
+}
+
+func (c *MultiScalarMulFoldedEdgeCasesTest) Define(api frontend.API) error {
+	cr, err := NewCurve(api)
+	if err != nil {
+		return err
+	}
+	ps := make([]*G1Affine, len(c.Points))
+	for i := range c.Points {
+		ps[i] = &c.Points[i]
+	}
+	ss := make([]*emulated.Element[ScalarField], len(c.Scalars))
+	for i := range c.Scalars {
+		ss[i] = &c.Scalars[i]
+	}
+	res, err := cr.MultiScalarMul(ps, ss, algopts.WithFoldingScalarMul(), algopts.WithUseSafe())
+	if err != nil {
+		return err
+	}
+	cr.AssertIsEqual(res, &c.Res)
+	return nil
+}
+
+func TestMultiScalarMulFoldedEdgeCases(t *testing.T) {
+	assert := test.NewAssert(t)
+	nbLen := 5
+	P := make([]bls12377.G1Affine, nbLen)
+	S := make([]fr.Element, nbLen)
+	S[0].SetOne()
+	S[1].SetRandom()
+	S[2].Square(&S[1])
+	S[3].Mul(&S[1], &S[2])
+	S[4].Mul(&S[1], &S[3])
+	for i := 0; i < nbLen; i++ {
+		P[i].ScalarMultiplicationBase(S[i].BigInt(new(big.Int)))
+	}
+	var res, infinity bls12377.G1Affine
+	_, err := res.MultiExp(P, S, ecc.MultiExpConfig{})
+
+	assert.NoError(err)
+	cP := make([]G1Affine, len(P))
+	cS := make([]emulated.Element[ScalarField], len(S))
+
+	// s^0 * (0,0) + s^1 * (0,0) + s^2 * (0,0) + s^3 * (0,0)  + s^4 * (0,0) == (0,0)
+	for i := range cP {
+		cP[i] = NewG1Affine(infinity)
+	}
+	// s0 = s
+	S[0].Set(&S[1])
+	for i := range cS {
+		cS[i] = NewScalar(S[i])
+	}
+	assignment1 := MultiScalarMulFoldedEdgeCasesTest{
+		Points:  cP,
+		Scalars: cS,
+		Res:     NewG1Affine(infinity),
+	}
+	err = test.IsSolved(&MultiScalarMulFoldedEdgeCasesTest{
+		Points:  make([]G1Affine, nbLen),
+		Scalars: make([]emulated.Element[ScalarField], nbLen),
+	}, &assignment1, ecc.BW6_761.ScalarField())
+	assert.NoError(err)
+
+	// 0^0 * P1 + 0 * P2 + 0 * P3 + 0 * P4 + 0 * P5 == P1
+	for i := range cP {
+		cP[i] = NewG1Affine(P[i])
+	}
+	for i := range cS {
+		cS[i] = emulated.ValueOf[emparams.BLS12377Fr](0)
+	}
+
+	assignment3 := MultiScalarMulFoldedEdgeCasesTest{
+		Points:  cP,
+		Scalars: cS,
+		Res:     NewG1Affine(P[0]),
+	}
+	err = test.IsSolved(&MultiScalarMulFoldedEdgeCasesTest{
+		Points:  make([]G1Affine, nbLen),
+		Scalars: make([]emulated.Element[ScalarField], nbLen),
+	}, &assignment3, ecc.BW6_761.ScalarField())
+	assert.NoError(err)
+}
+
+type MultiScalarMulFoldedTest struct {
+	Points  []G1Affine
+	Scalars []emulated.Element[ScalarField]
+	Res     G1Affine
+}
+
+func (c *MultiScalarMulFoldedTest) Define(api frontend.API) error {
+	cr, err := NewCurve(api)
+	if err != nil {
+		return err
+	}
+	ps := make([]*G1Affine, len(c.Points))
+	for i := range c.Points {
+		ps[i] = &c.Points[i]
+	}
+	ss := make([]*emulated.Element[ScalarField], len(c.Scalars))
+	for i := range c.Scalars {
+		ss[i] = &c.Scalars[i]
+	}
+	res, err := cr.MultiScalarMul(ps, ss, algopts.WithFoldingScalarMul())
+	if err != nil {
+		return err
+	}
+	cr.AssertIsEqual(res, &c.Res)
+	return nil
+}
+
+func TestMultiScalarMulFolded(t *testing.T) {
+	assert := test.NewAssert(t)
+	nbLen := 4
+	P := make([]bls12377.G1Affine, nbLen)
+	S := make([]fr.Element, nbLen)
+	// [s^0]P0 + [s^1]P1 + [s^2]P2 + [s^3]P3 = P0 + [s]P1 + [s^2]P2 + [s^3]P3
+	S[0].SetOne()
+	S[1].SetRandom()
+	S[2].Square(&S[1])
+	S[3].Mul(&S[1], &S[2])
+	for i := 0; i < nbLen; i++ {
+		P[i].ScalarMultiplicationBase(S[i].BigInt(new(big.Int)))
+	}
+	var res bls12377.G1Affine
+	_, err := res.MultiExp(P, S, ecc.MultiExpConfig{})
+
+	assert.NoError(err)
+	cP := make([]G1Affine, len(P))
+	for i := range cP {
+		cP[i] = NewG1Affine(P[i])
+	}
+	cS := make([]emulated.Element[ScalarField], len(S))
+	// s0 = s
+	S[0].Set(&S[1])
+	for i := range cS {
+		cS[i] = NewScalar(S[i])
+	}
+	assignment := MultiScalarMulFoldedTest{
+		Points:  cP,
+		Scalars: cS,
+		Res:     NewG1Affine(res),
+	}
+	err = test.IsSolved(&MultiScalarMulFoldedTest{
+		Points:  make([]G1Affine, nbLen),
+		Scalars: make([]emulated.Element[ScalarField], nbLen),
+	}, &assignment, ecc.BW6_761.ScalarField())
+	assert.NoError(err)
+}
diff --git a/std/algebra/native/sw_bls12377/g2.go b/std/algebra/native/sw_bls12377/g2.go
index 7d9e51689c..d183757a4d 100644
--- a/std/algebra/native/sw_bls12377/g2.go
+++ b/std/algebra/native/sw_bls12377/g2.go
@@ -24,14 +24,10 @@ import (
 
 	"github.com/consensys/gnark/constraint/solver"
 	"github.com/consensys/gnark/frontend"
+	"github.com/consensys/gnark/std/algebra/algopts"
 	"github.com/consensys/gnark/std/algebra/native/fields_bls12377"
 )
 
-// G2Jac point in Jacobian coords
-type G2Jac struct {
-	X, Y, Z fields_bls12377.E2
-}
-
 type g2AffP struct {
 	X, Y fields_bls12377.E2
 }
@@ -42,14 +38,6 @@ type G2Affine struct {
 	Lines *lineEvaluations
 }
 
-// Neg outputs -p
-func (p *G2Jac) Neg(api frontend.API, p1 G2Jac) *G2Jac {
-	p.Y.Neg(api, p1.Y)
-	p.X = p1.X
-	p.Z = p1.Z
-	return p
-}
-
 // Neg outputs -p
 func (p *g2AffP) Neg(api frontend.API, p1 g2AffP) *g2AffP {
 	p.Y.Neg(api, p1.Y)
@@ -83,86 +71,49 @@ func (p *g2AffP) AddAssign(api frontend.API, p1 g2AffP) *g2AffP {
 	return p
 }
 
-// AddAssign adds 2 point in Jacobian coordinates
-// p=p, a=p1
-func (p *G2Jac) AddAssign(api frontend.API, p1 *G2Jac) *G2Jac {
-
-	var Z1Z1, Z2Z2, U1, U2, S1, S2, H, I, J, r, V fields_bls12377.E2
-
-	Z1Z1.Square(api, p1.Z)
-
-	Z2Z2.Square(api, p.Z)
-
-	U1.Mul(api, p1.X, Z2Z2)
-
-	U2.Mul(api, p.X, Z1Z1)
-
-	S1.Mul(api, p1.Y, p.Z)
-	S1.Mul(api, S1, Z2Z2)
-
-	S2.Mul(api, p.Y, p1.Z)
-	S2.Mul(api, S2, Z1Z1)
-
-	H.Sub(api, U2, U1)
-
-	I.Add(api, H, H)
-	I.Square(api, I)
-
-	J.Mul(api, H, I)
-
-	r.Sub(api, S2, S1)
-	r.Add(api, r, r)
-
-	V.Mul(api, U1, I)
-
-	p.X.Square(api, r)
-	p.X.Sub(api, p.X, J)
-	p.X.Sub(api, p.X, V)
-	p.X.Sub(api, p.X, V)
-
-	p.Y.Sub(api, V, p.X)
-	p.Y.Mul(api, p.Y, r)
-
-	S1.Mul(api, J, S1)
-	S1.Add(api, S1, S1)
+func (p *g2AffP) AddUnified(api frontend.API, q g2AffP) *g2AffP {
+	// selector1 = 1 when p is (0,0) and 0 otherwise
+	selector1 := api.And(p.X.IsZero(api), p.Y.IsZero(api))
+	// selector2 = 1 when q is (0,0) and 0 otherwise
+	selector2 := api.And(q.X.IsZero(api), q.Y.IsZero(api))
+
+	// λ = ((p.x+q.x)² - p.x*q.x + a)/(p.y + q.y)
+	var pxqx, pxplusqx, num, denum, λ fields_bls12377.E2
+	pxqx.Mul(api, p.X, q.X)
+	pxplusqx.Add(api, p.X, q.X)
+	num.Mul(api, pxplusqx, pxplusqx)
+	num.Sub(api, num, pxqx)
+	denum.Add(api, p.Y, q.Y)
+	// if p.y + q.y = 0, assign dummy 1 to denum and continue
+	selector3 := denum.IsZero(api)
+	one := fields_bls12377.E2{A0: 1, A1: 0}
+	denum.Select(api, selector3, one, denum)
+	λ.DivUnchecked(api, num, denum)
+
+	// x = λ^2 - p.x - q.x
+	var xr, yr fields_bls12377.E2
+	xr.Square(api, λ)
+	xr.Sub(api, xr, pxplusqx)
+
+	// y = λ(p.x - xr) - p.y
+	yr.Sub(api, p.X, xr)
+	yr.Mul(api, yr, λ)
+	yr.Sub(api, yr, p.Y)
+	result := g2AffP{
+		X: xr,
+		Y: yr,
+	}
 
-	p.Y.Sub(api, p.Y, S1)
+	// if p=(0,0) return q
+	result.Select(api, selector1, q, result)
+	// if q=(0,0) return p
+	result.Select(api, selector2, *p, result)
+	// if p.y + q.y = 0, return (0, 0)
+	zero := fields_bls12377.E2{A0: 0, A1: 0}
+	result.Select(api, selector3, g2AffP{X: zero, Y: zero}, result)
 
-	p.Z.Add(api, p.Z, p1.Z)
-	p.Z.Square(api, p.Z)
-	p.Z.Sub(api, p.Z, Z1Z1)
-	p.Z.Sub(api, p.Z, Z2Z2)
-	p.Z.Mul(api, p.Z, H)
-
-	return p
-}
-
-// Double doubles a point in jacobian coords
-func (p *G2Jac) Double(api frontend.API, p1 G2Jac) *G2Jac {
-
-	var XX, YY, YYYY, ZZ, S, M, T fields_bls12377.E2
-
-	XX.Square(api, p.X)
-	YY.Square(api, p.Y)
-	YYYY.Square(api, YY)
-	ZZ.Square(api, p.Z)
-	S.Add(api, p.X, YY)
-	S.Square(api, S)
-	S.Sub(api, S, XX)
-	S.Sub(api, S, YYYY)
-	S.Add(api, S, S)
-	M.MulByFp(api, XX, 3) // M = 3*XX+a*ZZ², here a=0 (we suppose sw has j invariant 0)
-	p.Z.Add(api, p.Z, p.Y)
-	p.Z.Square(api, p.Z)
-	p.Z.Sub(api, p.Z, YY)
-	p.Z.Sub(api, p.Z, ZZ)
-	p.X.Square(api, M)
-	T.Add(api, S, S)
-	p.X.Sub(api, p.X, T)
-	p.Y.Sub(api, S, p.X)
-	p.Y.Mul(api, p.Y, M)
-	YYYY.MulByFp(api, YYYY, 8)
-	p.Y.Sub(api, p.Y, YYYY)
+	p.X = result.X
+	p.Y = result.Y
 
 	return p
 }
@@ -176,16 +127,6 @@ func (p *g2AffP) Select(api frontend.API, b frontend.Variable, p1, p2 g2AffP) *g
 	return p
 }
 
-// FromJac sets p to p1 in affine and returns it
-func (p *g2AffP) FromJac(api frontend.API, p1 G2Jac) *g2AffP {
-	var s fields_bls12377.E2
-	s.Mul(api, p1.Z, p1.Z)
-	p.X.DivUnchecked(api, p1.X, s)
-	s.Mul(api, s, p1.Z)
-	p.Y.DivUnchecked(api, p1.Y, s)
-	return p
-}
-
 // Double compute 2*p1, assign the result to p and return it
 // Only for curve with j invariant 0 (a=0).
 func (p *g2AffP) Double(api frontend.API, p1 g2AffP) *g2AffP {
@@ -219,11 +160,11 @@ func (p *g2AffP) Double(api frontend.API, p1 g2AffP) *g2AffP {
 // The method chooses an implementation based on scalar s. If it is constant,
 // then the compiled circuit depends on s. If it is variable type, then
 // the circuit is independent of the inputs.
-func (P *g2AffP) ScalarMul(api frontend.API, Q g2AffP, s interface{}) *g2AffP {
+func (P *g2AffP) ScalarMul(api frontend.API, Q g2AffP, s interface{}, opts ...algopts.AlgebraOption) *g2AffP {
 	if n, ok := api.Compiler().ConstantValue(s); ok {
-		return P.constScalarMul(api, Q, n)
+		return P.constScalarMul(api, Q, n, opts...)
 	} else {
-		return P.varScalarMul(api, Q, s)
+		return P.varScalarMul(api, Q, s, opts...)
 	}
 }
 
@@ -253,7 +194,11 @@ func init() {
 }
 
 // varScalarMul sets P = [s] Q and returns P.
-func (P *g2AffP) varScalarMul(api frontend.API, Q g2AffP, s frontend.Variable) *g2AffP {
+func (P *g2AffP) varScalarMul(api frontend.API, Q g2AffP, s frontend.Variable, opts ...algopts.AlgebraOption) *g2AffP {
+	cfg, err := algopts.NewConfig(opts...)
+	if err != nil {
+		panic(err)
+	}
 	// This method computes [s] Q. We use several methods to reduce the number
 	// of added constraints - first, instead of classical double-and-add, we use
 	// the optimized version from https://github.com/zcash/zcash/issues/3924
@@ -265,7 +210,13 @@ func (P *g2AffP) varScalarMul(api frontend.API, Q g2AffP, s frontend.Variable) *
 	// from a precomputed table. However, precomputing the table adds 12
 	// additional constraints and thus table-version is more expensive than
 	// addition-version.
-
+	var selector frontend.Variable
+	if cfg.UseSafe {
+		// if Q=(0,0) we assign a dummy (1,1) to Q and continue
+		selector = api.And(Q.X.IsZero(api), Q.Y.IsZero(api))
+		one := fields_bls12377.E2{A0: 1, A1: 0}
+		Q.Select(api, selector, g2AffP{X: one, Y: one}, Q)
+	}
 	// The context we are working is based on the `outer` curve. However, the
 	// points and the operations on the points are performed on the `inner`
 	// curve of the outer curve. We require some parameters from the inner
@@ -345,10 +296,22 @@ func (P *g2AffP) varScalarMul(api frontend.API, Q g2AffP, s frontend.Variable) *
 		Acc.DoubleAndAdd(api, &Acc, &B)
 	}
 
-	tableQ[0].AddAssign(api, Acc)
-	Acc.Select(api, s1bits[0], Acc, tableQ[0])
-	tablePhiQ[0].AddAssign(api, Acc)
-	Acc.Select(api, s2bits[0], Acc, tablePhiQ[0])
+	// i = 0
+	// When cfg.UseSafe is set, we use AddUnified instead of Add. This means
+	// when s=0 then Acc=(0,0) because AddUnified(Q, -Q) = (0,0).
+	if cfg.UseSafe {
+		tableQ[0].AddUnified(api, Acc)
+		Acc.Select(api, s1bits[0], Acc, tableQ[0])
+		tablePhiQ[0].AddUnified(api, Acc)
+		Acc.Select(api, s2bits[0], Acc, tablePhiQ[0])
+		zero := fields_bls12377.E2{A0: 0, A1: 0}
+		Acc.Select(api, selector, g2AffP{X: zero, Y: zero}, Acc)
+	} else {
+		tableQ[0].AddAssign(api, Acc)
+		Acc.Select(api, s1bits[0], Acc, tableQ[0])
+		tablePhiQ[0].AddAssign(api, Acc)
+		Acc.Select(api, s2bits[0], Acc, tablePhiQ[0])
+	}
 
 	P.X = Acc.X
 	P.Y = Acc.Y
@@ -357,7 +320,17 @@ func (P *g2AffP) varScalarMul(api frontend.API, Q g2AffP, s frontend.Variable) *
 }
 
 // constScalarMul sets P = [s] Q and returns P.
-func (P *g2AffP) constScalarMul(api frontend.API, Q g2AffP, s *big.Int) *g2AffP {
+func (P *g2AffP) constScalarMul(api frontend.API, Q g2AffP, s *big.Int, opts ...algopts.AlgebraOption) *g2AffP {
+	cfg, err := algopts.NewConfig(opts...)
+	if err != nil {
+		panic(err)
+	}
+	if s.BitLen() == 0 {
+		zero := fields_bls12377.E2{A0: 0, A1: 0}
+		P.X = zero
+		P.Y = zero
+		return P
+	}
 	// see the comments in varScalarMul. However, two-bit lookup is cheaper if
 	// bits are constant and here it makes sense to use the table in the main
 	// loop.
@@ -383,49 +356,60 @@ func (P *g2AffP) constScalarMul(api frontend.API, Q g2AffP, s *big.Int) *g2AffP
 	negPhiQ.Neg(api, phiQ)
 	var table [4]g2AffP
 	table[0] = negQ
-	table[0].AddAssign(api, negPhiQ)
 	table[1] = Q
-	table[1].AddAssign(api, negPhiQ)
 	table[2] = negQ
-	table[2].AddAssign(api, phiQ)
 	table[3] = Q
-	table[3].AddAssign(api, phiQ)
+
+	if cfg.UseSafe {
+		table[0].AddUnified(api, negPhiQ)
+		table[1].AddUnified(api, negPhiQ)
+		table[2].AddUnified(api, phiQ)
+		table[3].AddUnified(api, phiQ)
+	} else {
+		table[0].AddAssign(api, negPhiQ)
+		table[1].AddAssign(api, negPhiQ)
+		table[2].AddAssign(api, phiQ)
+		table[3].AddAssign(api, phiQ)
+	}
 
 	Acc = table[3]
 	// if both high bits are set, then we would get to the incomplete part,
 	// handle it separately.
 	if k[0].Bit(nbits-1) == 1 && k[1].Bit(nbits-1) == 1 {
-		Acc.Double(api, Acc)
-		Acc.AddAssign(api, table[3])
+		if cfg.UseSafe {
+			Acc.AddUnified(api, Acc)
+			Acc.AddUnified(api, table[3])
+		} else {
+			Acc.Double(api, Acc)
+			Acc.AddAssign(api, table[3])
+		}
 		nbits = nbits - 1
 	}
 	for i := nbits - 1; i > 0; i-- {
-		Acc.DoubleAndAdd(api, &Acc, &table[k[0].Bit(i)+2*k[1].Bit(i)])
+		if cfg.UseSafe {
+			Acc.AddUnified(api, Acc)
+			Acc.AddUnified(api, table[k[0].Bit(i)+2*k[1].Bit(i)])
+		} else {
+			Acc.DoubleAndAdd(api, &Acc, &table[k[0].Bit(i)+2*k[1].Bit(i)])
+		}
 	}
 
-	negQ.AddAssign(api, Acc)
-	Acc.Select(api, k[0].Bit(0), Acc, negQ)
-	negPhiQ.AddAssign(api, Acc)
+	// i = 0
+	if cfg.UseSafe {
+		negQ.AddUnified(api, Acc)
+		Acc.Select(api, k[0].Bit(0), Acc, negQ)
+		negPhiQ.AddUnified(api, Acc)
+	} else {
+		negQ.AddAssign(api, Acc)
+		Acc.Select(api, k[0].Bit(0), Acc, negQ)
+		negPhiQ.AddAssign(api, Acc)
+	}
 	Acc.Select(api, k[1].Bit(0), Acc, negPhiQ)
 	P.X, P.Y = Acc.X, Acc.Y
 
 	return P
 }
 
-// Assign a value to self (witness assignment)
-func (p *G2Jac) Assign(p1 *bls12377.G2Jac) {
-	p.X.Assign(&p1.X)
-	p.Y.Assign(&p1.Y)
-	p.Z.Assign(&p1.Z)
-}
-
-// AssertIsEqual constraint self to be equal to other into the given constraint system
-func (p *G2Jac) AssertIsEqual(api frontend.API, other G2Jac) {
-	p.X.AssertIsEqual(api, other.X)
-	p.Y.AssertIsEqual(api, other.Y)
-	p.Z.AssertIsEqual(api, other.Z)
-}
-
 // Assign a value to self (witness assignment)
 func (p *g2AffP) Assign(p1 *bls12377.G2Affine) {
 	p.X.Assign(&p1.X)
diff --git a/std/algebra/native/sw_bls12377/g2_test.go b/std/algebra/native/sw_bls12377/g2_test.go
index 7a6b941b63..c9ed38201b 100644
--- a/std/algebra/native/sw_bls12377/g2_test.go
+++ b/std/algebra/native/sw_bls12377/g2_test.go
@@ -23,48 +23,13 @@ import (
 	"github.com/consensys/gnark-crypto/ecc"
 	"github.com/consensys/gnark-crypto/ecc/bls12-377/fr"
 	"github.com/consensys/gnark/frontend"
+	"github.com/consensys/gnark/std/algebra/algopts"
+	"github.com/consensys/gnark/std/algebra/native/fields_bls12377"
 	"github.com/consensys/gnark/test"
 
 	bls12377 "github.com/consensys/gnark-crypto/ecc/bls12-377"
 )
 
-// -------------------------------------------------------------------------------------------------
-// Add jacobian
-
-type g2AddAssign struct {
-	A, B G2Jac
-	C    G2Jac `gnark:",public"`
-}
-
-func (circuit *g2AddAssign) Define(api frontend.API) error {
-	expected := circuit.A
-	expected.AddAssign(api, &circuit.B)
-	expected.AssertIsEqual(api, circuit.C)
-	return nil
-}
-
-func TestAddAssignG2(t *testing.T) {
-
-	// sample 2 random points
-	a := randomPointG2()
-	b := randomPointG2()
-
-	// create the cs
-	var circuit, witness g2AddAssign
-
-	// assign the inputs
-	witness.A.Assign(&a)
-	witness.B.Assign(&b)
-
-	// compute the result
-	a.AddAssign(&b)
-	witness.C.Assign(&a)
-
-	assert := test.NewAssert(t)
-	assert.CheckCircuit(&circuit, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_761))
-
-}
-
 // -------------------------------------------------------------------------------------------------
 // Add affine
 
@@ -106,41 +71,6 @@ func TestAddAssignAffineG2(t *testing.T) {
 
 }
 
-// -------------------------------------------------------------------------------------------------
-// Double Jacobian
-
-type g2DoubleAssign struct {
-	A G2Jac
-	C G2Jac `gnark:",public"`
-}
-
-func (circuit *g2DoubleAssign) Define(api frontend.API) error {
-	expected := circuit.A
-	expected.Double(api, circuit.A)
-	expected.AssertIsEqual(api, circuit.C)
-	return nil
-}
-
-func TestDoubleAssignG2(t *testing.T) {
-
-	// sample 2 random points
-	a := randomPointG2()
-
-	// create the cs
-	var circuit, witness g2DoubleAssign
-
-	// assign the inputs
-	witness.A.Assign(&a)
-
-	// compute the result
-	a.DoubleAssign()
-	witness.C.Assign(&a)
-
-	assert := test.NewAssert(t)
-	assert.CheckCircuit(&circuit, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_761))
-
-}
-
 // -------------------------------------------------------------------------------------------------
 // DoubleAndAdd affine
 
@@ -221,64 +151,72 @@ func TestDoubleAffineG2(t *testing.T) {
 }
 
 // -------------------------------------------------------------------------------------------------
-// Neg
+// Scalar multiplication
 
-type g2Neg struct {
-	A G2Jac
-	C G2Jac `gnark:",public"`
+type g2constantScalarMul struct {
+	A g2AffP
+	C g2AffP `gnark:",public"`
+	R *big.Int
 }
 
-func (circuit *g2Neg) Define(api frontend.API) error {
-	expected := G2Jac{}
-	expected.Neg(api, circuit.A)
+func (circuit *g2constantScalarMul) Define(api frontend.API) error {
+	expected := g2AffP{}
+	expected.constScalarMul(api, circuit.A, circuit.R)
 	expected.AssertIsEqual(api, circuit.C)
 	return nil
 }
 
-func TestNegG2(t *testing.T) {
-
-	// sample 2 random points
-	a := randomPointG2()
+func TestConstantScalarMulG2(t *testing.T) {
+	// sample random point
+	_a := randomPointG2()
+	var a, c bls12377.G2Affine
+	a.FromJacobian(&_a)
 
 	// create the cs
-	var circuit, witness g2Neg
-
+	var circuit, witness g2constantScalarMul
+	var r fr.Element
+	_, _ = r.SetRandom()
 	// assign the inputs
 	witness.A.Assign(&a)
-
 	// compute the result
-	a.Neg(&a)
-	witness.C.Assign(&a)
+	br := new(big.Int)
+	r.BigInt(br)
+	// br is a circuit parameter
+	circuit.R = br
+	_a.ScalarMultiplication(&_a, br)
+	c.FromJacobian(&_a)
+	witness.C.Assign(&c)
 
 	assert := test.NewAssert(t)
 	assert.CheckCircuit(&circuit, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_761))
 
 }
 
-// -------------------------------------------------------------------------------------------------
-// Scalar multiplication
-
-type g2constantScalarMul struct {
+type g2constantScalarMulEdgeCases struct {
 	A g2AffP
-	C g2AffP `gnark:",public"`
 	R *big.Int
 }
 
-func (circuit *g2constantScalarMul) Define(api frontend.API) error {
-	expected := g2AffP{}
-	expected.constScalarMul(api, circuit.A, circuit.R)
-	expected.AssertIsEqual(api, circuit.C)
+func (circuit *g2constantScalarMulEdgeCases) Define(api frontend.API) error {
+	expected1 := g2AffP{}
+	expected2 := g2AffP{}
+	zero := fields_bls12377.E2{A0: 0, A1: 0}
+	infinity := g2AffP{X: zero, Y: zero}
+	expected1.constScalarMul(api, circuit.A, big.NewInt(0))
+	expected2.constScalarMul(api, infinity, circuit.R, algopts.WithUseSafe())
+	expected1.AssertIsEqual(api, infinity)
+	expected2.AssertIsEqual(api, infinity)
 	return nil
 }
 
-func TestConstantScalarMulG2(t *testing.T) {
+func TestConstantScalarMulG2EdgeCases(t *testing.T) {
 	// sample random point
 	_a := randomPointG2()
-	var a, c bls12377.G2Affine
+	var a bls12377.G2Affine
 	a.FromJacobian(&_a)
 
 	// create the cs
-	var circuit, witness g2constantScalarMul
+	var circuit, witness g2constantScalarMulEdgeCases
 	var r fr.Element
 	_, _ = r.SetRandom()
 	// assign the inputs
@@ -288,9 +226,6 @@ func TestConstantScalarMulG2(t *testing.T) {
 	r.BigInt(br)
 	// br is a circuit parameter
 	circuit.R = br
-	_a.ScalarMultiplication(&_a, br)
-	c.FromJacobian(&_a)
-	witness.C.Assign(&c)
 
 	assert := test.NewAssert(t)
 	assert.CheckCircuit(&circuit, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_761))
@@ -333,6 +268,41 @@ func TestVarScalarMulG2(t *testing.T) {
 	assert.CheckCircuit(&circuit, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_761))
 }
 
+type g2varScalarMulEdgeCases struct {
+	A g2AffP
+	R frontend.Variable
+}
+
+func (circuit *g2varScalarMulEdgeCases) Define(api frontend.API) error {
+	expected1 := g2AffP{}
+	expected2 := g2AffP{}
+	zero := fields_bls12377.E2{A0: 0, A1: 0}
+	infinity := g2AffP{X: zero, Y: zero}
+	expected1.varScalarMul(api, circuit.A, 0, algopts.WithUseSafe())
+	expected2.varScalarMul(api, infinity, circuit.R, algopts.WithUseSafe())
+	expected1.AssertIsEqual(api, infinity)
+	expected2.AssertIsEqual(api, infinity)
+	return nil
+}
+
+func TestVarScalarMulG2EdgeCases(t *testing.T) {
+	// sample random point
+	_a := randomPointG2()
+	var a bls12377.G2Affine
+	a.FromJacobian(&_a)
+
+	// create the cs
+	var circuit, witness g2varScalarMulEdgeCases
+	var r fr.Element
+	_, _ = r.SetRandom()
+	witness.R = r.String()
+	// assign the inputs
+	witness.A.Assign(&a)
+
+	assert := test.NewAssert(t)
+	assert.CheckCircuit(&circuit, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_761))
+}
+
 type g2ScalarMul struct {
 	A    g2AffP
 	C    g2AffP `gnark:",public"`
diff --git a/std/algebra/native/sw_bls12377/pairing2.go b/std/algebra/native/sw_bls12377/pairing2.go
index 99cce9f5cf..e78bd8e738 100644
--- a/std/algebra/native/sw_bls12377/pairing2.go
+++ b/std/algebra/native/sw_bls12377/pairing2.go
@@ -107,7 +107,7 @@ func (c *Curve) jointScalarMul(P1, P2 *G1Affine, s1, s2 *Scalar, opts ...algopts
 	res := &G1Affine{}
 	varScalar1 := c.packScalarToVar(s1)
 	varScalar2 := c.packScalarToVar(s2)
-	res.jointScalarMul(c.api, *P1, *P2, varScalar1, varScalar2)
+	res.jointScalarMul(c.api, *P1, *P2, varScalar1, varScalar2, opts...)
 	return res
 }
 
@@ -119,7 +119,7 @@ func (c *Curve) ScalarMul(P *G1Affine, s *Scalar, opts ...algopts.AlgebraOption)
 		Y: P.Y,
 	}
 	varScalar := c.packScalarToVar(s)
-	res.ScalarMul(c.api, *P, varScalar)
+	res.ScalarMul(c.api, *P, varScalar, opts...)
 	return res
 }
 
@@ -128,7 +128,7 @@ func (c *Curve) ScalarMul(P *G1Affine, s *Scalar, opts ...algopts.AlgebraOption)
 func (c *Curve) ScalarMulBase(s *Scalar, opts ...algopts.AlgebraOption) *G1Affine {
 	res := new(G1Affine)
 	varScalar := c.packScalarToVar(s)
-	res.ScalarMulBase(c.api, varScalar)
+	res.ScalarMulBase(c.api, varScalar, opts...)
 	return res
 }
 
@@ -146,6 +146,10 @@ func (c *Curve) MultiScalarMul(P []*G1Affine, scalars []*Scalar, opts ...algopts
 	if err != nil {
 		return nil, fmt.Errorf("new config: %w", err)
 	}
+	addFn := c.Add
+	if cfg.UseSafe {
+		addFn = c.AddUnified
+	}
 	if !cfg.FoldMulti {
 		if len(P) != len(scalars) {
 			return nil, fmt.Errorf("mismatching points and scalars slice lengths")
@@ -160,7 +164,7 @@ func (c *Curve) MultiScalarMul(P []*G1Affine, scalars []*Scalar, opts ...algopts
 		}
 		for i := 1; i < n-1; i += 2 {
 			q := c.jointScalarMul(P[i-1], P[i], scalars[i-1], scalars[i], opts...)
-			res = c.Add(res, q)
+			res = addFn(res, q)
 		}
 		return res, nil
 	} else {
@@ -183,12 +187,12 @@ func (c *Curve) MultiScalarMul(P []*G1Affine, scalars []*Scalar, opts ...algopts
 
 		// points and scalars must be non-zero
 		var res G1Affine
-		res.scalarBitsMul(c.api, *P[len(P)-1], gamma1Bits, gamma2Bits)
+		res.scalarBitsMul(c.api, *P[len(P)-1], gamma1Bits, gamma2Bits, opts...)
 		for i := len(P) - 2; i > 0; i-- {
-			res = *c.Add(P[i], &res)
-			res.scalarBitsMul(c.api, res, gamma1Bits, gamma2Bits)
+			res = *addFn(P[i], &res)
+			res.scalarBitsMul(c.api, res, gamma1Bits, gamma2Bits, opts...)
 		}
-		res = *c.Add(P[0], &res)
+		res = *addFn(P[0], &res)
 		return &res, nil
 	}
 }
diff --git a/std/algebra/native/sw_bls24315/g1.go b/std/algebra/native/sw_bls24315/g1.go
index 8e06e311f8..3b958de711 100644
--- a/std/algebra/native/sw_bls24315/g1.go
+++ b/std/algebra/native/sw_bls24315/g1.go
@@ -27,24 +27,11 @@ import (
 	"github.com/consensys/gnark/std/algebra/algopts"
 )
 
-// G1Jac point in Jacobian coords
-type G1Jac struct {
-	X, Y, Z frontend.Variable
-}
-
 // G1Affine point in affine coords
 type G1Affine struct {
 	X, Y frontend.Variable
 }
 
-// Neg outputs -p
-func (p *G1Jac) Neg(api frontend.API, p1 G1Jac) *G1Jac {
-	p.X = p1.X
-	p.Y = api.Sub(0, p1.Y)
-	p.Z = p1.Z
-	return p
-}
-
 // Neg outputs -p
 func (p *G1Affine) Neg(api frontend.API, p1 G1Affine) *G1Affine {
 	p.X = p1.X
@@ -54,16 +41,16 @@ func (p *G1Affine) Neg(api frontend.API, p1 G1Affine) *G1Affine {
 
 func (p *G1Affine) AddUnified(api frontend.API, q G1Affine) *G1Affine {
 	// selector1 = 1 when p is (0,0) and 0 otherwise
-	selector1 := api.And(api.IsZero(&p.X), api.IsZero(&p.Y))
+	selector1 := api.And(api.IsZero(p.X), api.IsZero(p.Y))
 	// selector2 = 1 when q is (0,0) and 0 otherwise
-	selector2 := api.And(api.IsZero(&q.X), api.IsZero(&q.Y))
+	selector2 := api.And(api.IsZero(q.X), api.IsZero(q.Y))
 
 	// λ = ((p.x+q.x)² - p.x*q.x + a)/(p.y + q.y)
-	pxqx := api.Mul(&p.X, &q.X)
-	pxplusqx := api.Add(&p.X, &q.X)
+	pxqx := api.Mul(p.X, q.X)
+	pxplusqx := api.Add(p.X, q.X)
 	num := api.Mul(pxplusqx, pxplusqx)
 	num = api.Sub(num, pxqx)
-	denum := api.Add(&p.Y, &q.Y)
+	denum := api.Add(p.Y, q.Y)
 	// if p.y + q.y = 0, assign dummy 1 to denum and continue
 	selector3 := api.IsZero(denum)
 	denum = api.Select(selector3, 1, denum)
@@ -74,9 +61,9 @@ func (p *G1Affine) AddUnified(api frontend.API, q G1Affine) *G1Affine {
 	xr = api.Sub(xr, pxplusqx)
 
 	// y = λ(p.x - xr) - p.y
-	yr := api.Sub(&p.X, xr)
+	yr := api.Sub(p.X, xr)
 	yr = api.Mul(yr, λ)
-	yr = api.Sub(yr, &p.Y)
+	yr = api.Sub(yr, p.Y)
 	result := G1Affine{
 		X: xr,
 		Y: yr,
@@ -112,91 +99,6 @@ func (p *G1Affine) AddAssign(api frontend.API, p1 G1Affine) *G1Affine {
 	return p
 }
 
-// AddAssign adds 2 point in Jacobian coordinates
-// p=p, a=p1
-func (p *G1Jac) AddAssign(api frontend.API, p1 G1Jac) *G1Jac {
-
-	// get some Element from our pool
-	var Z1Z1, Z2Z2, U1, U2, S1, S2, H, I, J, r, V frontend.Variable
-
-	Z1Z1 = api.Mul(p1.Z, p1.Z)
-
-	Z2Z2 = api.Mul(p.Z, p.Z)
-
-	U1 = api.Mul(p1.X, Z2Z2)
-
-	U2 = api.Mul(p.X, Z1Z1)
-
-	S1 = api.Mul(p1.Y, p.Z)
-	S1 = api.Mul(S1, Z2Z2)
-
-	S2 = api.Mul(p.Y, p1.Z)
-	S2 = api.Mul(S2, Z1Z1)
-
-	H = api.Sub(U2, U1)
-
-	I = api.Add(H, H)
-	I = api.Mul(I, I)
-
-	J = api.Mul(H, I)
-
-	r = api.Sub(S2, S1)
-	r = api.Add(r, r)
-
-	V = api.Mul(U1, I)
-
-	p.X = api.Mul(r, r)
-	p.X = api.Sub(p.X, J)
-	p.X = api.Sub(p.X, V)
-	p.X = api.Sub(p.X, V)
-
-	p.Y = api.Sub(V, p.X)
-	p.Y = api.Mul(p.Y, r)
-
-	S1 = api.Mul(J, S1)
-	S1 = api.Add(S1, S1)
-
-	p.Y = api.Sub(p.Y, S1)
-
-	p.Z = api.Add(p.Z, p1.Z)
-	p.Z = api.Mul(p.Z, p.Z)
-	p.Z = api.Sub(p.Z, Z1Z1)
-	p.Z = api.Sub(p.Z, Z2Z2)
-	p.Z = api.Mul(p.Z, H)
-
-	return p
-}
-
-// DoubleAssign doubles the receiver point in jacobian coords and returns it
-func (p *G1Jac) DoubleAssign(api frontend.API) *G1Jac {
-	// get some Element from our pool
-	var XX, YY, YYYY, ZZ, S, M, T frontend.Variable
-
-	XX = api.Mul(p.X, p.X)
-	YY = api.Mul(p.Y, p.Y)
-	YYYY = api.Mul(YY, YY)
-	ZZ = api.Mul(p.Z, p.Z)
-	S = api.Add(p.X, YY)
-	S = api.Mul(S, S)
-	S = api.Sub(S, XX)
-	S = api.Sub(S, YYYY)
-	S = api.Add(S, S)
-	M = api.Mul(XX, 3) // M = 3*XX+a*ZZ², here a=0 (we suppose sw has j invariant 0)
-	p.Z = api.Add(p.Z, p.Y)
-	p.Z = api.Mul(p.Z, p.Z)
-	p.Z = api.Sub(p.Z, YY)
-	p.Z = api.Sub(p.Z, ZZ)
-	p.X = api.Mul(M, M)
-	T = api.Add(S, S)
-	p.X = api.Sub(p.X, T)
-	p.Y = api.Sub(S, p.X)
-	p.Y = api.Mul(p.Y, M)
-	YYYY = api.Mul(YYYY, 8)
-	p.Y = api.Sub(p.Y, YYYY)
-
-	return p
-}
-
 // Select sets p1 if b=1, p2 if b=0, and returns it. b must be boolean constrained
 func (p *G1Affine) Select(api frontend.API, b frontend.Variable, p1, p2 G1Affine) *G1Affine {
 
@@ -207,12 +109,19 @@ func (p *G1Affine) Select(api frontend.API, b frontend.Variable, p1, p2 G1Affine
 
 }
 
-// FromJac sets p to p1 in affine and returns it
-func (p *G1Affine) FromJac(api frontend.API, p1 G1Jac) *G1Affine {
-	s := api.Mul(p1.Z, p1.Z)
-	p.X = api.DivUnchecked(p1.X, s)
-	p.Y = api.DivUnchecked(p1.Y, api.Mul(s, p1.Z))
+// Lookup2 performs a 2-bit lookup between p1, p2, p3, p4 based on bits b0  and b1.
+// Returns:
+//   - p1 if b0=0 and b1=0,
+//   - p2 if b0=1 and b1=0,
+//   - p3 if b0=0 and b1=1,
+//   - p4 if b0=1 and b1=1.
+func (p *G1Affine) Lookup2(api frontend.API, b1, b2 frontend.Variable, p1, p2, p3, p4 G1Affine) *G1Affine {
+
+	p.X = api.Lookup2(b1, b2, p1.X, p2.X, p3.X, p4.X)
+	p.Y = api.Lookup2(b1, b2, p1.Y, p2.Y, p3.Y, p4.Y)
+
 	return p
+
 }
 
 // Double double a point in affine coords
@@ -242,11 +151,11 @@ func (p *G1Affine) Double(api frontend.API, p1 G1Affine) *G1Affine {
 // The method chooses an implementation based on scalar s. If it is constant,
 // then the compiled circuit depends on s. If it is variable type, then
 // the circuit is independent of the inputs.
-func (P *G1Affine) ScalarMul(api frontend.API, Q G1Affine, s interface{}) *G1Affine {
+func (P *G1Affine) ScalarMul(api frontend.API, Q G1Affine, s interface{}, opts ...algopts.AlgebraOption) *G1Affine {
 	if n, ok := api.Compiler().ConstantValue(s); ok {
-		return P.constScalarMul(api, Q, n)
+		return P.constScalarMul(api, Q, n, opts...)
 	} else {
-		return P.varScalarMul(api, Q, s)
+		return P.varScalarMul(api, Q, s, opts...)
 	}
 }
 
@@ -276,7 +185,11 @@ func init() {
 }
 
 // varScalarMul sets P = [s] Q and returns P.
-func (P *G1Affine) varScalarMul(api frontend.API, Q G1Affine, s frontend.Variable) *G1Affine {
+func (P *G1Affine) varScalarMul(api frontend.API, Q G1Affine, s frontend.Variable, opts ...algopts.AlgebraOption) *G1Affine {
+	cfg, err := algopts.NewConfig(opts...)
+	if err != nil {
+		panic(err)
+	}
 	// This method computes [s] Q. We use several methods to reduce the number
 	// of added constraints - first, instead of classical double-and-add, we use
 	// the optimized version from https://github.com/zcash/zcash/issues/3924
@@ -288,7 +201,12 @@ func (P *G1Affine) varScalarMul(api frontend.API, Q G1Affine, s frontend.Variabl
 	// from a precomputed table. However, precomputing the table adds 12
 	// additional constraints and thus table-version is more expensive than
 	// addition-version.
-
+	var selector frontend.Variable
+	if cfg.UseSafe {
+		// if Q=(0,0) we assign a dummy (1,1) to Q and continue
+		selector = api.And(api.IsZero(Q.X), api.IsZero(Q.Y))
+		Q.Select(api, selector, G1Affine{X: 1, Y: 1}, Q)
+	}
 	// The context we are working is based on the `outer` curve. However, the
 	// points and the operations on the points are performed on the `inner`
 	// curve of the outer curve. We require some parameters from the inner
@@ -368,10 +286,21 @@ func (P *G1Affine) varScalarMul(api frontend.API, Q G1Affine, s frontend.Variabl
 		Acc.DoubleAndAdd(api, &Acc, &B)
 	}
 
-	tableQ[0].AddAssign(api, Acc)
-	Acc.Select(api, s1bits[0], Acc, tableQ[0])
-	tablePhiQ[0].AddAssign(api, Acc)
-	Acc.Select(api, s2bits[0], Acc, tablePhiQ[0])
+	// i = 0
+	// When cfg.UseSafe is set, we use AddUnified instead of Add. This means
+	// when s=0 then Acc=(0,0) because AddUnified(Q, -Q) = (0,0).
+	if cfg.UseSafe {
+		tableQ[0].AddUnified(api, Acc)
+		Acc.Select(api, s1bits[0], Acc, tableQ[0])
+		tablePhiQ[0].AddUnified(api, Acc)
+		Acc.Select(api, s2bits[0], Acc, tablePhiQ[0])
+		Acc.Select(api, selector, G1Affine{X: 0, Y: 0}, Acc)
+	} else {
+		tableQ[0].AddAssign(api, Acc)
+		Acc.Select(api, s1bits[0], Acc, tableQ[0])
+		tablePhiQ[0].AddAssign(api, Acc)
+		Acc.Select(api, s2bits[0], Acc, tablePhiQ[0])
+	}
 
 	P.X = Acc.X
 	P.Y = Acc.Y
@@ -380,7 +309,16 @@ func (P *G1Affine) varScalarMul(api frontend.API, Q G1Affine, s frontend.Variabl
 }
 
 // constScalarMul sets P = [s] Q and returns P.
-func (P *G1Affine) constScalarMul(api frontend.API, Q G1Affine, s *big.Int) *G1Affine {
+func (P *G1Affine) constScalarMul(api frontend.API, Q G1Affine, s *big.Int, opts ...algopts.AlgebraOption) *G1Affine {
+	cfg, err := algopts.NewConfig(opts...)
+	if err != nil {
+		panic(err)
+	}
+	if s.BitLen() == 0 {
+		P.X = 0
+		P.Y = 0
+		return P
+	}
 	// see the comments in varScalarMul. However, two-bit lookup is cheaper if
 	// bits are constant and here it makes sense to use the table in the main
 	// loop.
@@ -406,49 +344,60 @@ func (P *G1Affine) constScalarMul(api frontend.API, Q G1Affine, s *big.Int) *G1A
 	negPhiQ.Neg(api, phiQ)
 	var table [4]G1Affine
 	table[0] = negQ
-	table[0].AddAssign(api, negPhiQ)
 	table[1] = Q
-	table[1].AddAssign(api, negPhiQ)
 	table[2] = negQ
-	table[2].AddAssign(api, phiQ)
 	table[3] = Q
-	table[3].AddAssign(api, phiQ)
+
+	if cfg.UseSafe {
+		table[0].AddUnified(api, negPhiQ)
+		table[1].AddUnified(api, negPhiQ)
+		table[2].AddUnified(api, phiQ)
+		table[3].AddUnified(api, phiQ)
+	} else {
+		table[0].AddAssign(api, negPhiQ)
+		table[1].AddAssign(api, negPhiQ)
+		table[2].AddAssign(api, phiQ)
+		table[3].AddAssign(api, phiQ)
+	}
 
 	Acc = table[3]
 	// if both high bits are set, then we would get to the incomplete part,
 	// handle it separately.
 	if k[0].Bit(nbits-1) == 1 && k[1].Bit(nbits-1) == 1 {
-		Acc.Double(api, Acc)
-		Acc.AddAssign(api, table[3])
+		if cfg.UseSafe {
+			Acc.AddUnified(api, Acc)
+			Acc.AddUnified(api, table[3])
+		} else {
+			Acc.Double(api, Acc)
+			Acc.AddAssign(api, table[3])
+		}
 		nbits = nbits - 1
 	}
 	for i := nbits - 1; i > 0; i-- {
-		Acc.DoubleAndAdd(api, &Acc, &table[k[0].Bit(i)+2*k[1].Bit(i)])
+		if cfg.UseSafe {
+			Acc.AddUnified(api, Acc)
+			Acc.AddUnified(api, table[k[0].Bit(i)+2*k[1].Bit(i)])
+		} else {
+			Acc.DoubleAndAdd(api, &Acc, &table[k[0].Bit(i)+2*k[1].Bit(i)])
+		}
 	}
 
-	negQ.AddAssign(api, Acc)
-	Acc.Select(api, k[0].Bit(0), Acc, negQ)
-	negPhiQ.AddAssign(api, Acc)
+	// i = 0
+	if cfg.UseSafe {
+		negQ.AddUnified(api, Acc)
+		Acc.Select(api, k[0].Bit(0), Acc, negQ)
+		negPhiQ.AddUnified(api, Acc)
+	} else {
+		negQ.AddAssign(api, Acc)
+		Acc.Select(api, k[0].Bit(0), Acc, negQ)
+		negPhiQ.AddAssign(api, Acc)
+	}
 	Acc.Select(api, k[1].Bit(0), Acc, negPhiQ)
 	P.X, P.Y = Acc.X, Acc.Y
 
 	return P
 }
 
-// Assign a value to self (witness assignment)
-func (p *G1Jac) Assign(p1 *bls24315.G1Jac) {
-	p.X = (fr.Element)(p1.X)
-	p.Y = (fr.Element)(p1.Y)
-	p.Z = (fr.Element)(p1.Z)
-}
-
-// AssertIsEqual constraint self to be equal to other into the given constraint system
-func (p *G1Jac) AssertIsEqual(api frontend.API, other G1Jac) {
-	api.AssertIsEqual(p.X, other.X)
-	api.AssertIsEqual(p.Y, other.Y)
-	api.AssertIsEqual(p.Z, other.Z)
-}
-
 // Assign a value to self (witness assignment)
 func (p *G1Affine) Assign(p1 *bls24315.G1Affine) {
 	p.X = (fr.Element)(p1.X)
@@ -473,10 +422,9 @@ func (p *G1Affine) DoubleAndAdd(api frontend.API, p1, p2 *G1Affine) *G1Affine {
 	x3 = api.Sub(x3, p2.X)
 
 	// omit y3 computation
-	// compute lambda2 = -lambda1-2*y1/(x3-x1)
+	// compute lambda2 = lambda1+2*y1/(x3-x1)
 	l2 := api.DivUnchecked(api.Add(p1.Y, p1.Y), api.Sub(x3, p1.X))
 	l2 = api.Add(l2, l1)
-	l2 = api.Neg(l2)
 
 	// compute x4 =lambda2**2-x1-x3
 	x4 := api.Mul(l2, l2)
@@ -484,7 +432,7 @@ func (p *G1Affine) DoubleAndAdd(api frontend.API, p1, p2 *G1Affine) *G1Affine {
 	x4 = api.Sub(x4, x3)
 
 	// compute y4 = lambda2*(x1 - x4)-y1
-	y4 := api.Sub(p1.X, x4)
+	y4 := api.Sub(x4, p1.X)
 	y4 = api.Mul(l2, y4)
 	y4 = api.Sub(y4, p1.Y)
 
@@ -495,36 +443,13 @@ func (p *G1Affine) DoubleAndAdd(api frontend.API, p1, p2 *G1Affine) *G1Affine {
 }
 
 // ScalarMulBase computes s * g1 and returns it, where g1 is the fixed generator. It doesn't modify s.
-func (P *G1Affine) ScalarMulBase(api frontend.API, s frontend.Variable) *G1Affine {
-
-	points := getCurvePoints()
-
-	sBits := api.ToBinary(s, 253)
-
-	var res, tmp G1Affine
-
-	// i = 1, 2
-	// gm[0] = 3g, gm[1] = 5g, gm[2] = 7g
-	res.X = api.Lookup2(sBits[1], sBits[2], points.G1x, points.G1m[0][0], points.G1m[1][0], points.G1m[2][0])
-	res.Y = api.Lookup2(sBits[1], sBits[2], points.G1y, points.G1m[0][1], points.G1m[1][1], points.G1m[2][1])
-
-	for i := 3; i < 253; i++ {
-		// gm[i] = [2^i]g
-		tmp.X = res.X
-		tmp.Y = res.Y
-		tmp.AddAssign(api, G1Affine{points.G1m[i][0], points.G1m[i][1]})
-		res.Select(api, sBits[i], tmp, res)
+func (P *G1Affine) ScalarMulBase(api frontend.API, s frontend.Variable, opts ...algopts.AlgebraOption) *G1Affine {
+	_, _, g1aff, _ := bls24315.Generators()
+	generator := G1Affine{
+		X: g1aff.X.BigInt(new(big.Int)),
+		Y: g1aff.Y.BigInt(new(big.Int)),
 	}
-
-	// i = 0
-	tmp.Neg(api, G1Affine{points.G1x, points.G1y})
-	tmp.AddAssign(api, res)
-	res.Select(api, sBits[0], res, tmp)
-
-	P.X = res.X
-	P.Y = res.Y
-
-	return P
+	return P.ScalarMul(api, generator, s, opts...)
 }
 
 // P = [s]Q + [t]R using Shamir's trick
@@ -613,8 +538,9 @@ func (P *G1Affine) jointScalarMul(api frontend.API, Q, R G1Affine, s, t frontend
 	return P
 }
 
-// scalarBitsMul...
-func (P *G1Affine) scalarBitsMul(api frontend.API, Q G1Affine, s1bits, s2bits []frontend.Variable, opts ...algopts.AlgebraOption) *G1Affine {
+// scalarBitsMul computes s * p and returns it where sBits is the bit decomposition of s. It doesn't modify p nor sBits.
+// The method is similar to varScalarMul.
+func (P *G1Affine) scalarBitsMul(api frontend.API, Q G1Affine, s1bits, s2bits []frontend.Variable) *G1Affine {
 	cc := getInnerCurveConfig(api.Compiler().Field())
 	nbits := cc.lambda.BitLen() + 1
 	var Acc /*accumulator*/, B, B2 /*tmp vars*/ G1Affine
diff --git a/std/algebra/native/sw_bls24315/g1_test.go b/std/algebra/native/sw_bls24315/g1_test.go
index f196347c89..1df72cdebe 100644
--- a/std/algebra/native/sw_bls24315/g1_test.go
+++ b/std/algebra/native/sw_bls24315/g1_test.go
@@ -24,6 +24,7 @@ import (
 	"github.com/consensys/gnark-crypto/ecc/bls24-315/fp"
 	"github.com/consensys/gnark-crypto/ecc/bls24-315/fr"
 	"github.com/consensys/gnark/frontend"
+	"github.com/consensys/gnark/std/algebra/algopts"
 	"github.com/consensys/gnark/std/math/emulated"
 	"github.com/consensys/gnark/test"
 
@@ -76,9 +77,7 @@ func (c *MarshalG1Test) Define(api frontend.API) error {
 	if err != nil {
 		return err
 	}
-	// we want to get the same output as gnark-crypto's marshal.
-	// It's a point on bls12-377 so the number of bytes is 96, as the
-	// field of definition of bls12-377 is 48 bytes long.
+	// the bits are layed out exactly as in gnark-crypto
 	r := ec.MarshalG1(c.P)
 	for i := range c.R {
 		api.AssertIsEqual(r[i], c.R[i])
@@ -118,43 +117,6 @@ func TestMarshalG1(t *testing.T) {
 	})
 }
 
-// -------------------------------------------------------------------------------------------------
-// Add jacobian
-
-type g1AddAssign struct {
-	A, B G1Jac
-	C    G1Jac `gnark:",public"`
-}
-
-func (circuit *g1AddAssign) Define(api frontend.API) error {
-	expected := circuit.A
-	expected.AddAssign(api, circuit.B)
-	expected.AssertIsEqual(api, circuit.C)
-	return nil
-}
-
-func TestAddAssignG1(t *testing.T) {
-
-	// sample 2 random points
-	a := randomPointG1()
-	b := randomPointG1()
-
-	// create the cs
-	var circuit, witness g1AddAssign
-
-	// assign the inputs
-	witness.A.Assign(&a)
-	witness.B.Assign(&b)
-
-	// compute the result
-	a.AddAssign(&b)
-	witness.C.Assign(&a)
-
-	assert := test.NewAssert(t)
-	assert.CheckCircuit(&circuit, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_633), test.NoProverChecks())
-
-}
-
 // -------------------------------------------------------------------------------------------------
 // Add affine
 
@@ -196,41 +158,6 @@ func TestAddAssignAffineG1(t *testing.T) {
 
 }
 
-// -------------------------------------------------------------------------------------------------
-// Double Jacobian
-
-type g1DoubleAssign struct {
-	A G1Jac
-	C G1Jac `gnark:",public"`
-}
-
-func (circuit *g1DoubleAssign) Define(api frontend.API) error {
-	expected := circuit.A
-	expected.DoubleAssign(api)
-	expected.AssertIsEqual(api, circuit.C)
-	return nil
-}
-
-func TestDoubleAssignG1(t *testing.T) {
-
-	// sample 2 random points
-	a := randomPointG1()
-
-	// create the cs
-	var circuit, witness g1DoubleAssign
-
-	// assign the inputs
-	witness.A.Assign(&a)
-
-	// compute the result
-	a.DoubleAssign()
-	witness.C.Assign(&a)
-
-	assert := test.NewAssert(t)
-	assert.CheckCircuit(&circuit, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_633), test.NoProverChecks())
-
-}
-
 // -------------------------------------------------------------------------------------------------
 // Double affine
 
@@ -307,37 +234,6 @@ func TestDoubleAndAddAffineG1(t *testing.T) {
 
 }
 
-// -------------------------------------------------------------------------------------------------
-// Neg
-
-type g1Neg struct {
-	A G1Jac
-	C G1Jac `gnark:",public"`
-}
-
-func (circuit *g1Neg) Define(api frontend.API) error {
-	expected := G1Jac{}
-	expected.Neg(api, circuit.A)
-	expected.AssertIsEqual(api, circuit.C)
-	return nil
-}
-
-func TestNegG1(t *testing.T) {
-
-	// sample 2 random points
-	a := randomPointG1()
-
-	// assign the inputs
-	var witness g1Neg
-	witness.A.Assign(&a)
-	a.Neg(&a)
-	witness.C.Assign(&a)
-
-	assert := test.NewAssert(t)
-	assert.CheckCircuit(&g1Neg{}, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_633), test.NoProverChecks())
-
-}
-
 // -------------------------------------------------------------------------------------------------
 // Scalar multiplication
 
@@ -380,6 +276,45 @@ func TestConstantScalarMulG1(t *testing.T) {
 
 }
 
+type g1constantScalarMulEdgeCases struct {
+	A G1Affine
+	R *big.Int
+}
+
+func (circuit *g1constantScalarMulEdgeCases) Define(api frontend.API) error {
+	expected1 := G1Affine{}
+	expected2 := G1Affine{}
+	infinity := G1Affine{X: 0, Y: 0}
+	expected1.constScalarMul(api, circuit.A, big.NewInt(0))
+	expected2.constScalarMul(api, infinity, circuit.R, algopts.WithUseSafe())
+	expected1.AssertIsEqual(api, infinity)
+	expected2.AssertIsEqual(api, infinity)
+	return nil
+}
+
+func TestConstantScalarMulG1EdgeCases(t *testing.T) {
+	// sample random point
+	_a := randomPointG1()
+	var a bls24315.G1Affine
+	a.FromJacobian(&_a)
+
+	// create the cs
+	var circuit, witness g1constantScalarMulEdgeCases
+	var r fr.Element
+	_, _ = r.SetRandom()
+	// assign the inputs
+	witness.A.Assign(&a)
+	// compute the result
+	br := new(big.Int)
+	r.BigInt(br)
+	// br is a circuit parameter
+	circuit.R = br
+
+	assert := test.NewAssert(t)
+	assert.CheckCircuit(&circuit, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_633))
+
+}
+
 type g1varScalarMul struct {
 	A G1Affine
 	C G1Affine `gnark:",public"`
@@ -456,6 +391,40 @@ func TestScalarMulG1(t *testing.T) {
 	assert.CheckCircuit(&circuit, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_633), test.NoProverChecks())
 }
 
+type g1varScalarMulEdgeCases struct {
+	A G1Affine
+	R frontend.Variable
+}
+
+func (circuit *g1varScalarMulEdgeCases) Define(api frontend.API) error {
+	expected1 := G1Affine{}
+	expected2 := G1Affine{}
+	infinity := G1Affine{X: 0, Y: 0}
+	expected1.varScalarMul(api, circuit.A, 0, algopts.WithUseSafe())
+	expected2.varScalarMul(api, infinity, circuit.R, algopts.WithUseSafe())
+	expected1.AssertIsEqual(api, infinity)
+	expected2.AssertIsEqual(api, infinity)
+	return nil
+}
+
+func TestVarScalarMulG1EdgeCases(t *testing.T) {
+	// sample random point
+	_a := randomPointG1()
+	var a bls24315.G1Affine
+	a.FromJacobian(&_a)
+
+	// create the cs
+	var circuit, witness g1varScalarMulEdgeCases
+	var r fr.Element
+	_, _ = r.SetRandom()
+	witness.R = r.String()
+	// assign the inputs
+	witness.A.Assign(&a)
+
+	assert := test.NewAssert(t)
+	assert.CheckCircuit(&circuit, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_633))
+}
+
 type g1varScalarMulBase struct {
 	C G1Affine `gnark:",public"`
 	R frontend.Variable
diff --git a/std/algebra/native/sw_bls24315/g2.go b/std/algebra/native/sw_bls24315/g2.go
index c7c23696df..70a5c31645 100644
--- a/std/algebra/native/sw_bls24315/g2.go
+++ b/std/algebra/native/sw_bls24315/g2.go
@@ -24,14 +24,10 @@ import (
 
 	"github.com/consensys/gnark/constraint/solver"
 	"github.com/consensys/gnark/frontend"
+	"github.com/consensys/gnark/std/algebra/algopts"
 	"github.com/consensys/gnark/std/algebra/native/fields_bls24315"
 )
 
-// G2Jac point in Jacobian coords
-type G2Jac struct {
-	X, Y, Z fields_bls24315.E4
-}
-
 type g2AffP struct {
 	X, Y fields_bls24315.E4
 }
@@ -42,14 +38,6 @@ type G2Affine struct {
 	Lines *lineEvaluations
 }
 
-// Neg outputs -p
-func (p *G2Jac) Neg(api frontend.API, p1 G2Jac) *G2Jac {
-	p.Y.Neg(api, p1.Y)
-	p.X = p1.X
-	p.Z = p1.Z
-	return p
-}
-
 // Neg outputs -p
 func (p *g2AffP) Neg(api frontend.API, p1 g2AffP) *g2AffP {
 	p.Y.Neg(api, p1.Y)
@@ -83,86 +71,49 @@ func (p *g2AffP) AddAssign(api frontend.API, p1 g2AffP) *g2AffP {
 	return p
 }
 
-// AddAssign adds 2 point in Jacobian coordinates
-// p=p, a=p1
-func (p *G2Jac) AddAssign(api frontend.API, p1 *G2Jac) *G2Jac {
-
-	var Z1Z1, Z2Z2, U1, U2, S1, S2, H, I, J, r, V fields_bls24315.E4
-
-	Z1Z1.Square(api, p1.Z)
-
-	Z2Z2.Square(api, p.Z)
-
-	U1.Mul(api, p1.X, Z2Z2)
-
-	U2.Mul(api, p.X, Z1Z1)
-
-	S1.Mul(api, p1.Y, p.Z)
-	S1.Mul(api, S1, Z2Z2)
-
-	S2.Mul(api, p.Y, p1.Z)
-	S2.Mul(api, S2, Z1Z1)
-
-	H.Sub(api, U2, U1)
-
-	I.Add(api, H, H)
-	I.Square(api, I)
-
-	J.Mul(api, H, I)
-
-	r.Sub(api, S2, S1)
-	r.Add(api, r, r)
-
-	V.Mul(api, U1, I)
-
-	p.X.Square(api, r)
-	p.X.Sub(api, p.X, J)
-	p.X.Sub(api, p.X, V)
-	p.X.Sub(api, p.X, V)
-
-	p.Y.Sub(api, V, p.X)
-	p.Y.Mul(api, p.Y, r)
-
-	S1.Mul(api, J, S1)
-	S1.Add(api, S1, S1)
-
-	p.Y.Sub(api, p.Y, S1)
-
-	p.Z.Add(api, p.Z, p1.Z)
-	p.Z.Square(api, p.Z)
-	p.Z.Sub(api, p.Z, Z1Z1)
-	p.Z.Sub(api, p.Z, Z2Z2)
-	p.Z.Mul(api, p.Z, H)
+func (p *g2AffP) AddUnified(api frontend.API, q g2AffP) *g2AffP {
+	// selector1 = 1 when p is (0,0) and 0 otherwise
+	selector1 := api.And(p.X.IsZero(api), p.Y.IsZero(api))
+	// selector2 = 1 when q is (0,0) and 0 otherwise
+	selector2 := api.And(q.X.IsZero(api), q.Y.IsZero(api))
+
+	// λ = ((p.x+q.x)² - p.x*q.x + a)/(p.y + q.y)
+	var pxqx, pxplusqx, num, denum, λ fields_bls24315.E4
+	pxqx.Mul(api, p.X, q.X)
+	pxplusqx.Add(api, p.X, q.X)
+	num.Mul(api, pxplusqx, pxplusqx)
+	num.Sub(api, num, pxqx)
+	denum.Add(api, p.Y, q.Y)
+	// if p.y + q.y = 0, assign dummy 1 to denum and continue
+	selector3 := denum.IsZero(api)
+	one := fields_bls24315.E4{B0: fields_bls24315.E2{A0: 1, A1: 0}, B1: fields_bls24315.E2{A0: 0, A1: 0}}
+	denum.Select(api, selector3, one, denum)
+	λ.DivUnchecked(api, num, denum)
+
+	// x = λ^2 - p.x - q.x
+	var xr, yr fields_bls24315.E4
+	xr.Square(api, λ)
+	xr.Sub(api, xr, pxplusqx)
+
+	// y = λ(p.x - xr) - p.y
+	yr.Sub(api, p.X, xr)
+	yr.Mul(api, yr, λ)
+	yr.Sub(api, yr, p.Y)
+	result := g2AffP{
+		X: xr,
+		Y: yr,
+	}
 
-	return p
-}
+	// if p=(0,0) return q
+	result.Select(api, selector1, q, result)
+	// if q=(0,0) return p
+	result.Select(api, selector2, *p, result)
+	// if p.y + q.y = 0, return (0, 0)
+	zero := fields_bls24315.E4{B0: fields_bls24315.E2{A0: 0, A1: 0}, B1: fields_bls24315.E2{A0: 0, A1: 0}}
+	result.Select(api, selector3, g2AffP{X: zero, Y: zero}, result)
 
-// Double doubles a point in jacobian coords
-func (p *G2Jac) Double(api frontend.API, p1 G2Jac) *G2Jac {
-
-	var XX, YY, YYYY, ZZ, S, M, T fields_bls24315.E4
-
-	XX.Square(api, p.X)
-	YY.Square(api, p.Y)
-	YYYY.Square(api, YY)
-	ZZ.Square(api, p.Z)
-	S.Add(api, p.X, YY)
-	S.Square(api, S)
-	S.Sub(api, S, XX)
-	S.Sub(api, S, YYYY)
-	S.Add(api, S, S)
-	M.MulByFp(api, XX, 3) // M = 3*XX+a*ZZ², here a=0 (we suppose sw has j invariant 0)
-	p.Z.Add(api, p.Z, p.Y)
-	p.Z.Square(api, p.Z)
-	p.Z.Sub(api, p.Z, YY)
-	p.Z.Sub(api, p.Z, ZZ)
-	p.X.Square(api, M)
-	T.Add(api, S, S)
-	p.X.Sub(api, p.X, T)
-	p.Y.Sub(api, S, p.X)
-	p.Y.Mul(api, p.Y, M)
-	YYYY.MulByFp(api, YYYY, 8)
-	p.Y.Sub(api, p.Y, YYYY)
+	p.X = result.X
+	p.Y = result.Y
 
 	return p
 }
@@ -176,16 +127,6 @@ func (p *g2AffP) Select(api frontend.API, b frontend.Variable, p1, p2 g2AffP) *g
 	return p
 }
 
-// FromJac sets p to p1 in affine and returns it
-func (p *g2AffP) FromJac(api frontend.API, p1 G2Jac) *g2AffP {
-	var s fields_bls24315.E4
-	s.Mul(api, p1.Z, p1.Z)
-	p.X.DivUnchecked(api, p1.X, s)
-	s.Mul(api, s, p1.Z)
-	p.Y.DivUnchecked(api, p1.Y, s)
-	return p
-}
-
 // Double compute 2*p1, assign the result to p and return it
 // Only for curve with j invariant 0 (a=0).
 func (p *g2AffP) Double(api frontend.API, p1 g2AffP) *g2AffP {
@@ -219,11 +160,11 @@ func (p *g2AffP) Double(api frontend.API, p1 g2AffP) *g2AffP {
 // The method chooses an implementation based on scalar s. If it is constant,
 // then the compiled circuit depends on s. If it is variable type, then
 // the circuit is independent of the inputs.
-func (P *g2AffP) ScalarMul(api frontend.API, Q g2AffP, s interface{}) *g2AffP {
+func (P *g2AffP) ScalarMul(api frontend.API, Q g2AffP, s interface{}, opts ...algopts.AlgebraOption) *g2AffP {
 	if n, ok := api.Compiler().ConstantValue(s); ok {
-		return P.constScalarMul(api, Q, n)
+		return P.constScalarMul(api, Q, n, opts...)
 	} else {
-		return P.varScalarMul(api, Q, s)
+		return P.varScalarMul(api, Q, s, opts...)
 	}
 }
 
@@ -253,7 +194,11 @@ func init() {
 }
 
 // varScalarMul sets P = [s] Q and returns P.
-func (P *g2AffP) varScalarMul(api frontend.API, Q g2AffP, s frontend.Variable) *g2AffP {
+func (P *g2AffP) varScalarMul(api frontend.API, Q g2AffP, s frontend.Variable, opts ...algopts.AlgebraOption) *g2AffP {
+	cfg, err := algopts.NewConfig(opts...)
+	if err != nil {
+		panic(err)
+	}
 	// This method computes [s] Q. We use several methods to reduce the number
 	// of added constraints - first, instead of classical double-and-add, we use
 	// the optimized version from https://github.com/zcash/zcash/issues/3924
@@ -265,7 +210,13 @@ func (P *g2AffP) varScalarMul(api frontend.API, Q g2AffP, s frontend.Variable) *
 	// from a precomputed table. However, precomputing the table adds 12
 	// additional constraints and thus table-version is more expensive than
 	// addition-version.
-
+	var selector frontend.Variable
+	if cfg.UseSafe {
+		// if Q=(0,0) we assign a dummy (1,1) to Q and continue
+		selector = api.And(Q.X.IsZero(api), Q.Y.IsZero(api))
+		one := fields_bls24315.E4{B0: fields_bls24315.E2{A0: 1, A1: 0}, B1: fields_bls24315.E2{A0: 0, A1: 0}}
+		Q.Select(api, selector, g2AffP{X: one, Y: one}, Q)
+	}
 	// The context we are working is based on the `outer` curve. However, the
 	// points and the operations on the points are performed on the `inner`
 	// curve of the outer curve. We require some parameters from the inner
@@ -345,10 +296,22 @@ func (P *g2AffP) varScalarMul(api frontend.API, Q g2AffP, s frontend.Variable) *
 		Acc.DoubleAndAdd(api, &Acc, &B)
 	}
 
-	tableQ[0].AddAssign(api, Acc)
-	Acc.Select(api, s1bits[0], Acc, tableQ[0])
-	tablePhiQ[0].AddAssign(api, Acc)
-	Acc.Select(api, s2bits[0], Acc, tablePhiQ[0])
+	// i = 0
+	// When cfg.UseSafe is set, we use AddUnified instead of Add. This means
+	// when s=0 then Acc=(0,0) because AddUnified(Q, -Q) = (0,0).
+	if cfg.UseSafe {
+		tableQ[0].AddUnified(api, Acc)
+		Acc.Select(api, s1bits[0], Acc, tableQ[0])
+		tablePhiQ[0].AddUnified(api, Acc)
+		Acc.Select(api, s2bits[0], Acc, tablePhiQ[0])
+		zero := fields_bls24315.E4{B0: fields_bls24315.E2{A0: 0, A1: 0}, B1: fields_bls24315.E2{A0: 0, A1: 0}}
+		Acc.Select(api, selector, g2AffP{X: zero, Y: zero}, Acc)
+	} else {
+		tableQ[0].AddAssign(api, Acc)
+		Acc.Select(api, s1bits[0], Acc, tableQ[0])
+		tablePhiQ[0].AddAssign(api, Acc)
+		Acc.Select(api, s2bits[0], Acc, tablePhiQ[0])
+	}
 
 	P.X = Acc.X
 	P.Y = Acc.Y
@@ -357,7 +320,17 @@ func (P *g2AffP) varScalarMul(api frontend.API, Q g2AffP, s frontend.Variable) *
 }
 
 // constScalarMul sets P = [s] Q and returns P.
-func (P *g2AffP) constScalarMul(api frontend.API, Q g2AffP, s *big.Int) *g2AffP {
+func (P *g2AffP) constScalarMul(api frontend.API, Q g2AffP, s *big.Int, opts ...algopts.AlgebraOption) *g2AffP {
+	cfg, err := algopts.NewConfig(opts...)
+	if err != nil {
+		panic(err)
+	}
+	if s.BitLen() == 0 {
+		zero := fields_bls24315.E4{B0: fields_bls24315.E2{A0: 0, A1: 0}, B1: fields_bls24315.E2{A0: 0, A1: 0}}
+		P.X = zero
+		P.Y = zero
+		return P
+	}
 	// see the comments in varScalarMul. However, two-bit lookup is cheaper if
 	// bits are constant and here it makes sense to use the table in the main
 	// loop.
@@ -366,7 +339,15 @@ func (P *g2AffP) constScalarMul(api frontend.API, Q g2AffP, s *big.Int) *g2AffP
 	s.Mod(s, cc.fr)
 	cc.phi2(api, &phiQ, &Q)
 
-	k := ecc.SplitScalar(s, cc.glvBasis)
+	var k [2]big.Int
+	// if s=0, assign dummy 1s to k[0] and k[1]
+	if s.BitLen() == 0 {
+		k[0].SetInt64(1)
+		k[1].SetInt64(1)
+	} else {
+		k = ecc.SplitScalar(s, cc.glvBasis)
+	}
+
 	if k[0].Sign() == -1 {
 		k[0].Neg(&k[0])
 		Q.Neg(api, Q)
@@ -383,49 +364,60 @@ func (P *g2AffP) constScalarMul(api frontend.API, Q g2AffP, s *big.Int) *g2AffP
 	negPhiQ.Neg(api, phiQ)
 	var table [4]g2AffP
 	table[0] = negQ
-	table[0].AddAssign(api, negPhiQ)
 	table[1] = Q
-	table[1].AddAssign(api, negPhiQ)
 	table[2] = negQ
-	table[2].AddAssign(api, phiQ)
 	table[3] = Q
-	table[3].AddAssign(api, phiQ)
+
+	if cfg.UseSafe {
+		table[0].AddUnified(api, negPhiQ)
+		table[1].AddUnified(api, negPhiQ)
+		table[2].AddUnified(api, phiQ)
+		table[3].AddUnified(api, phiQ)
+	} else {
+		table[0].AddAssign(api, negPhiQ)
+		table[1].AddAssign(api, negPhiQ)
+		table[2].AddAssign(api, phiQ)
+		table[3].AddAssign(api, phiQ)
+	}
 
 	Acc = table[3]
 	// if both high bits are set, then we would get to the incomplete part,
 	// handle it separately.
 	if k[0].Bit(nbits-1) == 1 && k[1].Bit(nbits-1) == 1 {
-		Acc.Double(api, Acc)
-		Acc.AddAssign(api, table[3])
+		if cfg.UseSafe {
+			Acc.AddUnified(api, Acc)
+			Acc.AddUnified(api, table[3])
+		} else {
+			Acc.Double(api, Acc)
+			Acc.AddAssign(api, table[3])
+		}
 		nbits = nbits - 1
 	}
 	for i := nbits - 1; i > 0; i-- {
-		Acc.DoubleAndAdd(api, &Acc, &table[k[0].Bit(i)+2*k[1].Bit(i)])
+		if cfg.UseSafe {
+			Acc.AddUnified(api, Acc)
+			Acc.AddUnified(api, table[k[0].Bit(i)+2*k[1].Bit(i)])
+		} else {
+			Acc.DoubleAndAdd(api, &Acc, &table[k[0].Bit(i)+2*k[1].Bit(i)])
+		}
 	}
 
-	negQ.AddAssign(api, Acc)
-	Acc.Select(api, k[0].Bit(0), Acc, negQ)
-	negPhiQ.AddAssign(api, Acc)
+	// i = 0
+	if cfg.UseSafe {
+		negQ.AddUnified(api, Acc)
+		Acc.Select(api, k[0].Bit(0), Acc, negQ)
+		negPhiQ.AddUnified(api, Acc)
+	} else {
+		negQ.AddAssign(api, Acc)
+		Acc.Select(api, k[0].Bit(0), Acc, negQ)
+		negPhiQ.AddAssign(api, Acc)
+	}
 	Acc.Select(api, k[1].Bit(0), Acc, negPhiQ)
 	P.X, P.Y = Acc.X, Acc.Y
 
 	return P
 }
 
-// Assign a value to self (witness assignment)
-func (p *G2Jac) Assign(p1 *bls24315.G2Jac) {
-	p.X.Assign(&p1.X)
-	p.Y.Assign(&p1.Y)
-	p.Z.Assign(&p1.Z)
-}
-
-// AssertIsEqual constraint self to be equal to other into the given constraint system
-func (p *G2Jac) AssertIsEqual(api frontend.API, other G2Jac) {
-	p.X.AssertIsEqual(api, other.X)
-	p.Y.AssertIsEqual(api, other.Y)
-	p.Z.AssertIsEqual(api, other.Z)
-}
-
 // Assign a value to self (witness assignment)
 func (p *g2AffP) Assign(p1 *bls24315.G2Affine) {
 	p.X.Assign(&p1.X)
diff --git a/std/algebra/native/sw_bls24315/g2_test.go b/std/algebra/native/sw_bls24315/g2_test.go
index 8f5ddfbfe1..879a454bae 100644
--- a/std/algebra/native/sw_bls24315/g2_test.go
+++ b/std/algebra/native/sw_bls24315/g2_test.go
@@ -23,48 +23,13 @@ import (
 	"github.com/consensys/gnark-crypto/ecc"
 	"github.com/consensys/gnark-crypto/ecc/bls24-315/fr"
 	"github.com/consensys/gnark/frontend"
+	"github.com/consensys/gnark/std/algebra/algopts"
+	"github.com/consensys/gnark/std/algebra/native/fields_bls24315"
 	"github.com/consensys/gnark/test"
 
 	bls24315 "github.com/consensys/gnark-crypto/ecc/bls24-315"
 )
 
-// -------------------------------------------------------------------------------------------------
-// Add jacobian
-
-type g2AddAssign struct {
-	A, B G2Jac
-	C    G2Jac `gnark:",public"`
-}
-
-func (circuit *g2AddAssign) Define(api frontend.API) error {
-	expected := circuit.A
-	expected.AddAssign(api, &circuit.B)
-	expected.AssertIsEqual(api, circuit.C)
-	return nil
-}
-
-func TestAddAssignG2(t *testing.T) {
-
-	// sample 2 random points
-	a := randomPointG2()
-	b := randomPointG2()
-
-	// create the cs
-	var circuit, witness g2AddAssign
-
-	// assign the inputs
-	witness.A.Assign(&a)
-	witness.B.Assign(&b)
-
-	// compute the result
-	a.AddAssign(&b)
-	witness.C.Assign(&a)
-
-	assert := test.NewAssert(t)
-	assert.CheckCircuit(&circuit, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_633), test.NoProverChecks())
-
-}
-
 // -------------------------------------------------------------------------------------------------
 // Add affine
 
@@ -106,41 +71,6 @@ func TestAddAssignAffineG2(t *testing.T) {
 
 }
 
-// -------------------------------------------------------------------------------------------------
-// Double Jacobian
-
-type g2DoubleAssign struct {
-	A G2Jac
-	C G2Jac `gnark:",public"`
-}
-
-func (circuit *g2DoubleAssign) Define(api frontend.API) error {
-	expected := circuit.A
-	expected.Double(api, circuit.A)
-	expected.AssertIsEqual(api, circuit.C)
-	return nil
-}
-
-func TestDoubleAssignG2(t *testing.T) {
-
-	// sample 2 random points
-	a := randomPointG2()
-
-	// create the cs
-	var circuit, witness g2DoubleAssign
-
-	// assign the inputs
-	witness.A.Assign(&a)
-
-	// compute the result
-	a.DoubleAssign()
-	witness.C.Assign(&a)
-
-	assert := test.NewAssert(t)
-	assert.CheckCircuit(&circuit, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_633), test.NoProverChecks())
-
-}
-
 // -------------------------------------------------------------------------------------------------
 // DoubleAndAdd affine
 
@@ -221,64 +151,72 @@ func TestDoubleAffineG2(t *testing.T) {
 }
 
 // -------------------------------------------------------------------------------------------------
-// Neg
+// Scalar multiplication
 
-type g2Neg struct {
-	A G2Jac
-	C G2Jac `gnark:",public"`
+type g2constantScalarMul struct {
+	A g2AffP
+	C g2AffP `gnark:",public"`
+	R *big.Int
 }
 
-func (circuit *g2Neg) Define(api frontend.API) error {
-	expected := G2Jac{}
-	expected.Neg(api, circuit.A)
+func (circuit *g2constantScalarMul) Define(api frontend.API) error {
+	expected := g2AffP{}
+	expected.constScalarMul(api, circuit.A, circuit.R)
 	expected.AssertIsEqual(api, circuit.C)
 	return nil
 }
 
-func TestNegG2(t *testing.T) {
-
-	// sample 2 random points
-	a := randomPointG2()
+func TestConstantScalarMulG2(t *testing.T) {
+	// sample random point
+	_a := randomPointG2()
+	var a, c bls24315.G2Affine
+	a.FromJacobian(&_a)
 
 	// create the cs
-	var circuit, witness g2Neg
-
+	var circuit, witness g2constantScalarMul
+	var r fr.Element
+	_, _ = r.SetRandom()
 	// assign the inputs
 	witness.A.Assign(&a)
-
 	// compute the result
-	a.Neg(&a)
-	witness.C.Assign(&a)
+	br := new(big.Int)
+	r.BigInt(br)
+	// br is a circuit parameter
+	circuit.R = br
+	_a.ScalarMultiplication(&_a, br)
+	c.FromJacobian(&_a)
+	witness.C.Assign(&c)
 
 	assert := test.NewAssert(t)
 	assert.CheckCircuit(&circuit, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_633), test.NoProverChecks())
 
 }
 
-// -------------------------------------------------------------------------------------------------
-// Scalar multiplication
-
-type g2constantScalarMul struct {
+type g2constantScalarMulEdgeCases struct {
 	A g2AffP
-	C g2AffP `gnark:",public"`
 	R *big.Int
 }
 
-func (circuit *g2constantScalarMul) Define(api frontend.API) error {
-	expected := g2AffP{}
-	expected.constScalarMul(api, circuit.A, circuit.R)
-	expected.AssertIsEqual(api, circuit.C)
+func (circuit *g2constantScalarMulEdgeCases) Define(api frontend.API) error {
+	expected1 := g2AffP{}
+	expected2 := g2AffP{}
+	zero := fields_bls24315.E4{B0: fields_bls24315.E2{A0: 0, A1: 0}, B1: fields_bls24315.E2{A0: 0, A1: 0}}
+	infinity := g2AffP{X: zero, Y: zero}
+	expected1.constScalarMul(api, circuit.A, big.NewInt(0))
+	expected2.constScalarMul(api, infinity, circuit.R, algopts.WithUseSafe())
+	expected1.AssertIsEqual(api, infinity)
+	expected2.AssertIsEqual(api, infinity)
 	return nil
 }
 
-func TestConstantScalarMulG2(t *testing.T) {
+func TestConstantScalarMulG2EdgeCases(t *testing.T) {
 	// sample random point
 	_a := randomPointG2()
-	var a, c bls24315.G2Affine
+	var a bls24315.G2Affine
 	a.FromJacobian(&_a)
 
 	// create the cs
-	var circuit, witness g2constantScalarMul
+	var circuit, witness g2constantScalarMulEdgeCases
 	var r fr.Element
 	_, _ = r.SetRandom()
 	// assign the inputs
@@ -288,13 +226,45 @@ func TestConstantScalarMulG2(t *testing.T) {
 	r.BigInt(br)
 	// br is a circuit parameter
 	circuit.R = br
-	_a.ScalarMultiplication(&_a, br)
-	c.FromJacobian(&_a)
-	witness.C.Assign(&c)
 
 	assert := test.NewAssert(t)
-	assert.CheckCircuit(&circuit, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_633), test.NoProverChecks())
+	assert.CheckCircuit(&circuit, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_633))
+
+}
+
+type g2varScalarMulEdgeCases struct {
+	A g2AffP
+	R frontend.Variable
+}
+
+func (circuit *g2varScalarMulEdgeCases) Define(api frontend.API) error {
+	expected1 := g2AffP{}
+	expected2 := g2AffP{}
+	zero := fields_bls24315.E4{B0: fields_bls24315.E2{A0: 0, A1: 0}, B1: fields_bls24315.E2{A0: 0, A1: 0}}
+	infinity := g2AffP{X: zero, Y: zero}
+	expected1.varScalarMul(api, circuit.A, 0, algopts.WithUseSafe())
+	expected2.varScalarMul(api, infinity, circuit.R, algopts.WithUseSafe())
+	expected1.AssertIsEqual(api, infinity)
+	expected2.AssertIsEqual(api, infinity)
+	return nil
+}
+
+func TestVarScalarMulG2EdgeCases(t *testing.T) {
+	// sample random point
+	_a := randomPointG2()
+	var a bls24315.G2Affine
+	a.FromJacobian(&_a)
+
+	// create the cs
+	var circuit, witness g2varScalarMulEdgeCases
+	var r fr.Element
+	_, _ = r.SetRandom()
+	witness.R = r.String()
+	// assign the inputs
+	witness.A.Assign(&a)
 
+	assert := test.NewAssert(t)
+	assert.CheckCircuit(&circuit, test.WithValidAssignment(&witness), test.WithCurves(ecc.BW6_633))
 }
 
 type g2varScalarMul struct {
diff --git a/std/recursion/plonk/verifier.go b/std/recursion/plonk/verifier.go
index 67c7764095..bcd83f13ba 100644
--- a/std/recursion/plonk/verifier.go
+++ b/std/recursion/plonk/verifier.go
@@ -1002,11 +1002,11 @@ func (v *Verifier[FR, G1El, G2El, GtEl]) PrepareVerification(vk VerifyingKey[FR,
 		_s1, _s2, // second & third part
 	)
 
-	linearizedPolynomialDigest, err := v.curve.MultiScalarMul(points, scalars)
+	linearizedPolynomialDigest, err := v.curve.MultiScalarMul(points, scalars /*, algopts.WithUseSafe()*/) // in PLONK Wo Commit ==> use algopts.WithUseSafe()
 	if err != nil {
 		return nil, nil, nil, fmt.Errorf("linearized polynomial digest MSM: %w", err)
 	}
-	linearizedPolynomialDigest = v.curve.Add(linearizedPolynomialDigest, &vk.Qk.G1El) // Qk=0 in PLONK W\ Commit ==> use AddUnified
+	linearizedPolynomialDigest = v.curve. /*AddUnified*/ Add(linearizedPolynomialDigest, &vk.Qk.G1El) // in PLONK Wo Commit ==> use AddUnified
 
 	// Fold the first proof
 	digestsToFold := make([]kzg.Commitment[G1El], len(vk.Qcp)+7)
diff --git a/std/recursion/plonk/verifier_test.go b/std/recursion/plonk/verifier_test.go
index f621a83f38..92a8016b54 100644
--- a/std/recursion/plonk/verifier_test.go
+++ b/std/recursion/plonk/verifier_test.go
@@ -42,10 +42,12 @@ func (c *OuterCircuit[FR, G1El, G2El, GtEl]) Define(api frontend.API) error {
 }
 
 /*
-TODO: Tests without api.Commit fail because the current optimized MultiScalarMul (JointScalarMul in particular) requires input points to be distinct.
-
-//-----------------------------------------------------------------
+///-----------------------------------------------------------------
 // Without api.Commit
+//
+// For this, in verifier.go, use:
+// - linearizedPolynomialDigest, err := v.curve.MultiScalarMul(points, scalars, algopts.WithUseSafe())
+// - linearizedPolynomialDigest = v.curve.AddUnified(linearizedPolynomialDigest, &vk.Qk.G1El)
 
 type InnerCircuitNativeWoCommit struct {
 	P, Q frontend.Variable
@@ -58,13 +60,13 @@ func (c *InnerCircuitNativeWoCommit) Define(api frontend.API) error {
 	return nil
 }
 
-func getInnerWoCommit(assert *test.Assert, field, outer *big.Int) (constraint.ConstraintSystem, plonk.VerifyingKey, witness.Witness, plonk.Proof) {
+func getInnerWoCommit(assert *test.Assert, field, outer *big.Int) (constraint.ConstraintSystem, native_plonk.VerifyingKey, witness.Witness, native_plonk.Proof) {
 	innerCcs, err := frontend.Compile(field, scs.NewBuilder, &InnerCircuitNativeWoCommit{})
 	assert.NoError(err)
 	srs, srsLagrange, err := unsafekzg.NewSRS(innerCcs)
 	assert.NoError(err)
 
-	innerPK, innerVK, err := plonk.Setup(innerCcs, srs, srsLagrange)
+	innerPK, innerVK, err := native_plonk.Setup(innerCcs, srs, srsLagrange)
 	assert.NoError(err)
 
 	// inner proof
@@ -75,11 +77,11 @@ func getInnerWoCommit(assert *test.Assert, field, outer *big.Int) (constraint.Co
 	}
 	innerWitness, err := frontend.NewWitness(innerAssignment, field)
 	assert.NoError(err)
-	innerProof, err := plonk.Prove(innerCcs, innerPK, innerWitness, GetNativeProverOptions(outer, field))
+	innerProof, err := native_plonk.Prove(innerCcs, innerPK, innerWitness, GetNativeProverOptions(outer, field))
 	assert.NoError(err)
 	innerPubWitness, err := innerWitness.Public()
 	assert.NoError(err)
-	err = plonk.Verify(innerProof, innerVK, innerPubWitness, GetNativeVerifierOptions(outer, field))
+	err = native_plonk.Verify(innerProof, innerVK, innerPubWitness, GetNativeVerifierOptions(outer, field))
 	assert.NoError(err)
 	return innerCcs, innerVK, innerPubWitness, innerProof
 }
@@ -100,12 +102,11 @@ func TestBLS12InBW6WoCommit(t *testing.T) {
 	outerCircuit := &OuterCircuit[sw_bls12377.ScalarField, sw_bls12377.G1Affine, sw_bls12377.G2Affine, sw_bls12377.GT]{
 		InnerWitness: PlaceholderWitness[sw_bls12377.ScalarField](innerCcs),
 		Proof:        PlaceholderProof[sw_bls12377.ScalarField, sw_bls12377.G1Affine, sw_bls12377.G2Affine](innerCcs),
-		VerifyingKey: PlaceholderVerifyingKey[sw_bls12377.ScalarField, sw_bls12377.G1Affine, sw_bls12377.G2Affine](innerCcs),
+		VerifyingKey: circuitVk,
 	}
 	outerAssignment := &OuterCircuit[sw_bls12377.ScalarField, sw_bls12377.G1Affine, sw_bls12377.G2Affine, sw_bls12377.GT]{
 		InnerWitness: circuitWitness,
 		Proof:        circuitProof,
-		VerifyingKey: circuitVk,
 	}
 	err = test.IsSolved(outerCircuit, outerAssignment, ecc.BW6_761.ScalarField())
 	assert.NoError(err)
@@ -127,12 +128,11 @@ func TestBW6InBN254WoCommit(t *testing.T) {
 	outerCircuit := &OuterCircuit[sw_bw6761.ScalarField, sw_bw6761.G1Affine, sw_bw6761.G2Affine, sw_bw6761.GTEl]{
 		InnerWitness: PlaceholderWitness[sw_bw6761.ScalarField](innerCcs),
 		Proof:        PlaceholderProof[sw_bw6761.ScalarField, sw_bw6761.G1Affine, sw_bw6761.G2Affine](innerCcs),
-		VerifyingKey: PlaceholderVerifyingKey[sw_bw6761.ScalarField, sw_bw6761.G1Affine, sw_bw6761.G2Affine](innerCcs),
+		VerifyingKey: circuitVk,
 	}
 	outerAssignment := &OuterCircuit[sw_bw6761.ScalarField, sw_bw6761.G1Affine, sw_bw6761.G2Affine, sw_bw6761.GTEl]{
 		InnerWitness: circuitWitness,
 		Proof:        circuitProof,
-		VerifyingKey: circuitVk,
 	}
 	err = test.IsSolved(outerCircuit, outerAssignment, ecc.BN254.ScalarField())
 	assert.NoError(err)
@@ -154,12 +154,11 @@ func TestBLS12381InBN254WoCommit(t *testing.T) {
 	outerCircuit := &OuterCircuit[sw_bls12381.ScalarField, sw_bls12381.G1Affine, sw_bls12381.G2Affine, sw_bls12381.GTEl]{
 		InnerWitness: PlaceholderWitness[sw_bls12381.ScalarField](innerCcs),
 		Proof:        PlaceholderProof[sw_bls12381.ScalarField, sw_bls12381.G1Affine, sw_bls12381.G2Affine](innerCcs),
-		VerifyingKey: PlaceholderVerifyingKey[sw_bls12381.ScalarField, sw_bls12381.G1Affine, sw_bls12381.G2Affine](innerCcs),
+		VerifyingKey: circuitVk,
 	}
 	outerAssignment := &OuterCircuit[sw_bls12381.ScalarField, sw_bls12381.G1Affine, sw_bls12381.G2Affine, sw_bls12381.GTEl]{
 		InnerWitness: circuitWitness,
 		Proof:        circuitProof,
-		VerifyingKey: circuitVk,
 	}
 	err = test.IsSolved(outerCircuit, outerAssignment, ecc.BN254.ScalarField())
 	assert.NoError(err)