diff --git a/packages/baby-jubjub/src/baby-jubjub.ts b/packages/baby-jubjub/src/baby-jubjub.ts
index 7f523d7e8..e6304b4d2 100644
--- a/packages/baby-jubjub/src/baby-jubjub.ts
+++ b/packages/baby-jubjub/src/baby-jubjub.ts
@@ -26,6 +26,8 @@ export const Base8: Point<bigint> = [
 const a = Fr.e(BigInt("168700"))
 const d = Fr.e(BigInt("168696"))
 
+export const id: Point<bigint> = [0n, 1n]
+
 // The Baby JubJub curve 'E(F_r)' is equal to the subgroup of 'F_r'-rational points of 'E'.
 export const order = BigInt("21888242871839275222246405745257275088614511777268538073601725287587578984328")
 export const subOrder = scalar.shiftRight(order, BigInt(3))
@@ -65,25 +67,38 @@ export function addPoint(p1: Point<bigint>, p2: Point<bigint>): Point<bigint> {
 /**
  * Performs a scalar multiplication by starting from the 'base' point and 'adding'
  * it to itself 'e' times.
+ * This algorithm is called 'Montgomery Ladder'. See {@link https://en.wikipedia.org/wiki/Elliptic_curve_point_multiplication#Montgomery_ladder}
+ * This works given the following invariant: At each step, R0 will be r_0*base where r_0 is the prefix of e
+ * written in binary and R1 will be (r_0+1)*base. In other words: at iteration i of the loop, r_0's binary
+ * representation will be the first i+1 most significant bits of e. If the upcoming bit is a 0, we just have to
+ * double R0 and add R0 to R1 to maintain the invariant. If it is a 1, we have to double R0 and add 1*base
+ * (or add R1, which is the same as (r_0+1)*base), and double R1 to maintain the invariant.
  * @param base The base point used as a starting point.
  * @param e A secret number representing the private key.
  * @returns The resulting point representing the public key.
  */
 export function mulPointEscalar(base: Point<bigint>, e: bigint): Point<bigint> {
-    let res: Point<bigint> = [Fr.e(BigInt(0)), Fr.e(BigInt(1))]
-    let rem: bigint = e
-    let exp: Point<bigint> = base
-
-    while (!scalar.isZero(rem)) {
-        if (scalar.isOdd(rem)) {
-            res = addPoint(res, exp)
+    e %= order
+    // set a bit above the maximum value so that the exponent
+    // variable will always be 254 bits for subsequent operations
+    // the 254th bit should be ignored in any operations below
+    e += 1n << 254n
+
+    let R0: Point<bigint> = id
+    let R1: Point<bigint> = base
+
+    // 'order' is a number of 254 bits, such as 1n<<253n. Therefore, we initialize the mask as 1<<253
+    for (let mask = 1n << 253n; mask > 0; mask >>= 1n) {
+        if (e & mask) {
+            R0 = addPoint(R0, R1)
+            R1 = addPoint(R1, R1)
+        } else {
+            R1 = addPoint(R0, R1)
+            R0 = addPoint(R0, R0)
         }
-
-        exp = addPoint(exp, exp)
-        rem = scalar.shiftRight(rem, BigInt(1))
     }
 
-    return res
+    return R0
 }
 
 /**
diff --git a/packages/baby-jubjub/tests/index.test.ts b/packages/baby-jubjub/tests/index.test.ts
index b1de333b9..ab51707d7 100644
--- a/packages/baby-jubjub/tests/index.test.ts
+++ b/packages/baby-jubjub/tests/index.test.ts
@@ -1,6 +1,19 @@
 import { babyjub } from "circomlibjs"
 import { utils } from "ffjavascript"
-import { Base8, Point, addPoint, inCurve, mulPointEscalar, packPoint, r, unpackPoint } from "../src"
+import * as scalar from "@zk-kit/utils/scalar"
+import {
+    Base8,
+    Point,
+    addPoint,
+    inCurve,
+    mulPointEscalar,
+    packPoint,
+    r,
+    unpackPoint,
+    order,
+    subOrder,
+    id
+} from "../src"
 import { tonelliShanks } from "../src/sqrt"
 
 describe("BabyJubjub", () => {
@@ -8,6 +21,63 @@ describe("BabyJubjub", () => {
 
     let publicKey: Point<bigint>
 
+    it("Test point addition and inCurve", async () => {
+        expect(inCurve(id)).toBeTruthy()
+        expect(inCurve([BigInt(1), BigInt(0)])).toBeFalsy()
+        expect(addPoint(id, id).toString()).toBe(id.toString())
+
+        const p1: Point<bigint> = [
+            BigInt("17777552123799933955779906779655732241715742912184938656739573121738514868268"),
+            BigInt("2626589144620713026669568689430873010625803728049924121243784502389097019475")
+        ]
+        const p2: Point<bigint> = [
+            BigInt("16540640123574156134436876038791482806971768689494387082833631921987005038935"),
+            BigInt("20819045374670962167435360035096875258406992893633759881276124905556507972311")
+        ]
+        const p3: Point<bigint> = [
+            BigInt("7916061937171219682591368294088513039687205273691143098332585753343424131937"),
+            BigInt("14035240266687799601661095864649209771790948434046947201833777492504781204499")
+        ]
+
+        expect(inCurve(p1)).toBeTruthy()
+        expect(inCurve(p2)).toBeTruthy()
+        expect(inCurve(p3)).toBeTruthy()
+
+        expect(addPoint(p1, p2).toString()).toBe(p3.toString())
+    })
+
+    it("Test point multiplication with small values", async () => {
+        const P: Point<bigint> = [
+            BigInt("17777552123799933955779906779655732241715742912184938656739573121738514868268"),
+            BigInt("2626589144620713026669568689430873010625803728049924121243784502389097019475")
+        ]
+
+        expect(mulPointEscalar(P, BigInt(0)).toString()).toBe(id.toString())
+        expect(mulPointEscalar(P, BigInt(1)).toString()).toBe(P.toString())
+        expect(mulPointEscalar(P, BigInt(2)).toString()).toBe(addPoint(P, P).toString())
+        expect(mulPointEscalar(P, BigInt(3)).toString()).toBe(addPoint(addPoint(P, P), P).toString())
+
+        expect(mulPointEscalar(id, BigInt(1)).toString()).toBe(id.toString())
+        expect(mulPointEscalar(id, BigInt(14134324)).toString()).toBe(id.toString())
+    })
+
+    it("Test base point order", async () => {
+        expect(scalar.shiftRight(order, BigInt(3))).toBe(subOrder)
+        const G: Point<bigint> = [
+            BigInt("995203441582195749578291179787384436505546430278305826713579947235728471134"),
+            BigInt("5472060717959818805561601436314318772137091100104008585924551046643952123905")
+        ]
+        const p1: Point<bigint> = mulPointEscalar(G, BigInt(8) * subOrder)
+        expect(p1.toString()).toBe(id.toString())
+        const p2 = mulPointEscalar(Base8, subOrder)
+        expect(p2.toString()).toBe(id.toString())
+
+        const random = BigInt("38275423985628165")
+        expect(mulPointEscalar(Base8, random).toString()).toBe(
+            mulPointEscalar(Base8, random + BigInt(543523) * subOrder).toString()
+        )
+    })
+
     it("Should add 1 point to the curve", async () => {
         const p1: Point<bigint> = [BigInt(0), BigInt(1)]