diff --git a/internal/edwards25519/scalarMul.go b/internal/edwards25519/scalarMul.go
index 050ebae..e1f4a85 100644
--- a/internal/edwards25519/scalarMul.go
+++ b/internal/edwards25519/scalarMul.go
@@ -113,7 +113,52 @@ func (v *ProjP3) VartimeDoubleBaseMul(a, b *scalar.Scalar, A *ProjP3) *ProjP3 {
 //
 // The multiscalar multiplication is performed in constant time.
 func (v *ProjP3) MultiscalarMul(scalars []scalar.Scalar, points []*ProjP3) *ProjP3 {
-	panic("unimplemented")
+	if len(scalars) != len(points) {
+		panic("called MultiscalarMul with different size inputs")
+	}
+
+	// Proceed as in the single-base case, but share doublings
+	// between each point in the multiscalar equation.
+
+	// Build lookup tables for each point
+	tables := make([]ProjLookupTable, len(points))
+	for i := range tables {
+		tables[i].FromP3(points[i])
+	}
+	// Compute signed radix-16 digits for each scalar
+	digits := make([][64]int8, len(scalars))
+	for i := range digits {
+		digits[i] = scalars[i].SignedRadix16()
+	}
+
+	// Unwrap first loop iteration to save computing 16*identity
+	multiple := &ProjCached{}
+	tmp1 := &ProjP1xP1{}
+	tmp2 := &ProjP2{}
+	// Lookup-and-add the appropriate multiple of each input point
+	for j := range tables {
+		tables[j].SelectInto(multiple, digits[j][63])
+		tmp1.Add(v, multiple) // tmp1 = v + x_(j,63)*Q in P1xP1 coords
+		v.FromP1xP1(tmp1)     // update v
+	}
+	tmp2.FromP3(v) // set up tmp2 = v in P2 coords for next iteration
+	for i := 62; i >= 0; i-- {
+		tmp1.Double(tmp2)    // tmp1 =  2*(prev) in P1xP1 coords
+		tmp2.FromP1xP1(tmp1) // tmp2 =  2*(prev) in P2 coords
+		tmp1.Double(tmp2)    // tmp1 =  4*(prev) in P1xP1 coords
+		tmp2.FromP1xP1(tmp1) // tmp2 =  4*(prev) in P2 coords
+		tmp1.Double(tmp2)    // tmp1 =  8*(prev) in P1xP1 coords
+		tmp2.FromP1xP1(tmp1) // tmp2 =  8*(prev) in P2 coords
+		tmp1.Double(tmp2)    // tmp1 = 16*(prev) in P1xP1 coords
+		v.FromP1xP1(tmp1)    //    v = 16*(prev) in P3 coords
+		// Lookup-and-add the appropriate multiple of each input point
+		for j := range tables {
+			tables[j].SelectInto(multiple, digits[j][i])
+			tmp1.Add(v, multiple) // tmp1 = v + x_(j,i)*Q in P1xP1 coords
+			v.FromP1xP1(tmp1)     // update v
+		}
+		tmp2.FromP3(v) // set up tmp2 = v in P2 coords for next iteration
+	}
 	return v
 }
 
diff --git a/internal/edwards25519/scalarMul_test.go b/internal/edwards25519/scalarMul_test.go
index 7a93437..989dd08 100644
--- a/internal/edwards25519/scalarMul_test.go
+++ b/internal/edwards25519/scalarMul_test.go
@@ -122,3 +122,27 @@ func TestScalarMulMatchesBasepointMul(t *testing.T) {
 		t.Error(err)
 	}
 }
+
+func TestMultiScalarMulMatchesBasepointMul(t *testing.T) {
+	multiScalarMulMatchesBasepointMul := func(x, y, z scalar.Scalar) bool {
+		// FIXME opaque scalars
+		x[31] &= 127
+		y[31] &= 127
+		z[31] &= 127
+		var p, q1, q2, q3, check ProjP3
+
+		p.MultiscalarMul([]scalar.Scalar{x, y, z}, []*ProjP3{&B, &B, &B})
+
+		q1.BasepointMul(&x)
+		q2.BasepointMul(&y)
+		q3.BasepointMul(&z)
+		check.Zero()
+		check.Add(&q1, &q2).Add(&check, &q3)
+
+		return p.Equal(&check) == 1
+	}
+
+	if err := quick.Check(multiScalarMulMatchesBasepointMul, quickCheckConfig); err != nil {
+		t.Error(err)
+	}
+}