diff --git a/bitmapcontainer.go b/bitmapcontainer.go
index b5b9d226..ee5f5fc3 100644
--- a/bitmapcontainer.go
+++ b/bitmapcontainer.go
@@ -798,8 +798,26 @@ func (bc *bitmapContainer) iand(a container) container {
 }
 
 func (bc *bitmapContainer) iandRun16(rc *runContainer16) container {
-	rcb := newBitmapContainerFromRun(rc)
-	return bc.iandBitmap(rcb)
+	if len(rc.iv) > runAndScratchIntervals {
+		// Clearing gap by gap costs a call per interval; past this many the
+		// run is built in scratch and ANDed in one flat pass instead.
+		var scratch [bitmapContainerSize]uint64
+		for i := range rc.iv {
+			setBitmapRange(scratch[:], int(rc.iv[i].start), int(rc.iv[i].last())+1)
+		}
+		return bc.iandBitmap(&bitmapContainer{bitmap: scratch[:]})
+	}
+	card := rc.andBitmapContainerCardinality(bc)
+	if card <= arrayDefaultMaxSize {
+		answer := newArrayContainerCapacity(card)
+		for i := range rc.iv {
+			answer.content = appendBitmapRange(answer.content, bc.bitmap, int(rc.iv[i].start), int(rc.iv[i].last())+1)
+		}
+		return answer
+	}
+	clearBitmapGaps(bc.bitmap, rc.iv, 0, maxCapacity)
+	bc.cardinality = card
+	return bc
 }
 
 func (bc *bitmapContainer) iandArray(ac *arrayContainer) container {
@@ -809,15 +827,13 @@ func (bc *bitmapContainer) iandArray(ac *arrayContainer) container {
 
 func (bc *bitmapContainer) andArray(value2 *arrayContainer) *arrayContainer {
 	answer := newArrayContainerCapacity(len(value2.content))
-	answer.content = answer.content[:cap(answer.content)]
-	c := value2.getCardinality()
+	out, bitmap := answer.content[:cap(answer.content)], bc.bitmap
 	pos := 0
-	for k := 0; k < c; k++ {
-		v := value2.content[k]
-		answer.content[pos] = v
-		pos += int(bc.bitValue(v))
+	for _, v := range value2.content {
+		out[pos] = v
+		pos += int((bitmap[v>>6] >> (v & 63)) & 1)
 	}
-	answer.content = answer.content[:pos]
+	answer.content = out[:pos]
 	return answer
 }
 
@@ -831,6 +847,26 @@ func (bc *bitmapContainer) andArrayCardinality(value2 *arrayContainer) int {
 	return pos
 }
 
+func (bc *bitmapContainer) intersectsRange(start, end uint) bool {
+	if start >= end {
+		return false
+	}
+	firstword, endword := start/64, (end-1)/64
+	lo, hi := ^uint64(0)<<(start%64), ^uint64(0)>>((64-end)&63)
+	if firstword == endword {
+		return bc.bitmap[firstword]&lo&hi != 0
+	}
+	if bc.bitmap[firstword]&lo != 0 || bc.bitmap[endword]&hi != 0 {
+		return true
+	}
+	for _, w := range bc.bitmap[firstword+1 : endword] {
+		if w != 0 {
+			return true
+		}
+	}
+	return false
+}
+
 func (bc *bitmapContainer) getCardinalityInRange(start, end uint) int {
 	if start >= end {
 		return 0
diff --git a/fastaggregation.go b/fastaggregation.go
index 3881b88a..6dfd2a16 100644
--- a/fastaggregation.go
+++ b/fastaggregation.go
@@ -125,26 +125,191 @@ func (x1 *Bitmap) repairAfterLazy() {
 	}
 }
 
-// FastAnd computes the intersection between many bitmaps quickly
-// Compared to the And function, it can take many bitmaps as input, thus saving the trouble
-// of manually calling "And" many times.
+// FastAnd computes the intersection between many bitmaps quickly.
+// Compared to the And function, it can take many bitmaps as input.
 //
-// Performance hints: if you have very large and tiny bitmaps,
-// it may be beneficial performance-wise to put a tiny bitmap
-// in first position.
+// With three or more inputs, the keys of the smallest input are walked
+// against every other input and probed with a container-level intersects
+// test; only a key every input shares a value on is intersected, and a key
+// of bitmaps is counted before its result is allocated.
 func FastAnd(bitmaps ...*Bitmap) *Bitmap {
-	if len(bitmaps) == 0 {
+	switch len(bitmaps) {
+	case 0:
 		return NewBitmap()
-	} else if len(bitmaps) == 1 {
+	case 1:
 		return bitmaps[0].Clone()
+	case 2:
+		return And(bitmaps[0], bitmaps[1])
 	}
-	answer := And(bitmaps[0], bitmaps[1])
-	for _, bm := range bitmaps[2:] {
-		answer.And(bm)
+	driver := 0
+	for i := 1; i < len(bitmaps); i++ {
+		if bitmaps[i].highlowcontainer.size() < bitmaps[driver].highlowcontainer.size() {
+			driver = i
+		}
+	}
+	dra := &bitmaps[driver].highlowcontainer
+	// A cursor per input and the containers of one key, on the stack for up
+	// to eight inputs; the containers are only made once a key needs them.
+	var posBuf [8]int
+	var csBuf [8]container
+	var pos []int
+	var cs []container
+	if len(bitmaps) <= len(posBuf) {
+		pos, cs = posBuf[:len(bitmaps)], csBuf[:len(bitmaps)]
+	} else {
+		pos = make([]int, len(bitmaps))
+	}
+	answer := NewBitmap()
+	// Keys come in order, so every input is walked, not searched, and a key
+	// is probed for a shared value with each input before its containers are
+	// intersected: a key with nothing in common allocates nothing.
+keys:
+	for j := 0; j < dra.size(); j++ {
+		key := dra.getKeyAtIndex(j)
+		dc := dra.getContainerAtIndex(j)
+		pos[driver] = j
+		for i, bm := range bitmaps {
+			if i == driver {
+				continue
+			}
+			ra := &bm.highlowcontainer
+			if pos[i] < ra.size() && ra.getKeyAtIndex(pos[i]) < key {
+				pos[i] = ra.advanceUntil(key, pos[i]) // searches from pos+1
+			}
+			if pos[i] >= ra.size() {
+				break keys
+			}
+			if ra.getKeyAtIndex(pos[i]) != key || !dc.intersects(ra.getContainerAtIndex(pos[i])) {
+				continue keys
+			}
+		}
+		if cs == nil {
+			cs = make([]container, len(bitmaps))
+		}
+		for i, bm := range bitmaps {
+			cs[i] = bm.highlowcontainer.getContainerAtIndex(pos[i])
+			pos[i]++
+		}
+		if c := andK(cs); c != nil {
+			answer.highlowcontainer.appendContainer(key, c, false)
+		}
 	}
 	return answer
 }
 
+// andK intersects the containers of one key, three or more, and returns nil
+// when the intersection is empty. cs is scratch and is reordered in place.
+func andK(cs []container) container {
+	// A full run is the identity and drops out; the smallest array leads
+	// and the rest keep the caller's order.
+	n, smallest := 0, -1
+	for _, c := range cs {
+		switch x := c.(type) {
+		case *runContainer16:
+			if x.isFull() {
+				continue
+			}
+		case *arrayContainer:
+			if smallest < 0 || x.getCardinality() < cs[smallest].getCardinality() {
+				smallest = n
+			}
+		}
+		cs[n] = c
+		n++
+	}
+	switch cs = cs[:n]; len(cs) {
+	case 0:
+		return newRunContainer16Range(0, maxCapacity-1)
+	case 1:
+		return cs[0].clone()
+	}
+	if smallest >= 0 {
+		a := cs[smallest]
+		copy(cs[1:smallest+1], cs[:smallest])
+		cs[0] = a
+		return andKChain(cs)
+	}
+	n = 0
+	for i, c := range cs {
+		if _, ok := c.(*runContainer16); ok {
+			cs[n], cs[i] = cs[i], cs[n]
+			n++
+		}
+	}
+	if n == len(cs) {
+		return andKChain(cs)
+	}
+	return andKBitmaps(cs[n:], cs[:n])
+}
+
+// andKChain intersects cs in order with the library's own kernels, the first
+// pair into a fresh container and the rest in place, so nothing bigger than
+// that first result is built.
+func andKChain(cs []container) container {
+	c := cs[0].and(cs[1])
+	for i := 2; i < len(cs) && !c.isEmpty(); i++ {
+		c = c.iand(cs[i])
+	}
+	if c.isEmpty() {
+		return nil
+	}
+	return c
+}
+
+// andKBitmaps ANDs the bitmaps into stack scratch, stopping at the first
+// empty prefix, over the words the runs leave: every run narrows the span to
+// its extent, and only a run with gaps is folded into a mask afterwards, so
+// a range built with AddRange costs no intersection at all.
+func andKBitmaps(bms, runs []container) container {
+	lo, hi := 0, maxCapacity-1
+	n := 0
+	for _, c := range runs {
+		rc := c.(*runContainer16)
+		lo, hi = max(lo, int(rc.minimum())), min(hi, int(rc.maximum()))
+		if len(rc.iv) > 1 {
+			runs[n] = c
+			n++
+		}
+	}
+	if lo > hi {
+		return nil
+	}
+	runs = runs[:n]
+	first, last := lo/64, hi/64
+	var scratch [bitmapContainerSize]uint64
+	w := scratch[first : last+1]
+	src := bms[0].(*bitmapContainer).bitmap[first : last+1]
+	card := uint64(0)
+	if len(bms) == 1 {
+		if card = andCardSlice(w, src, src); card == 0 {
+			return nil
+		}
+	}
+	for _, c := range bms[1:] {
+		if card = andCardSlice(w, src, c.(*bitmapContainer).bitmap[first:last+1]); card == 0 {
+			return nil
+		}
+		src = w
+	}
+	if len(runs) > 0 || lo%64 != 0 || (hi+1)%64 != 0 {
+		resetBitmapRange(scratch[:], first*64, lo)
+		resetBitmapRange(scratch[:], hi+1, (last+1)*64)
+		if len(runs) > 0 {
+			mask := runs[0].(*runContainer16)
+			for _, c := range runs[1:] {
+				if mask = mask.intersect(c.(*runContainer16)); len(mask.iv) == 0 {
+					return nil
+				}
+			}
+			clearBitmapGaps(scratch[:], mask.iv, lo, hi+1)
+		}
+		if card = popcntSlice(w); card == 0 {
+			return nil
+		}
+	}
+	return containerFromWords(scratch[:], first, last, int(card))
+}
+
 // FastOr computes the union between many bitmaps quickly, as opposed to having to call Or repeatedly.
 // It might also be faster than calling Or repeatedly.
 func FastOr(bitmaps ...*Bitmap) *Bitmap {
diff --git a/fastand_kway_test.go b/fastand_kway_test.go
new file mode 100644
index 00000000..f7102a33
--- /dev/null
+++ b/fastand_kway_test.go
@@ -0,0 +1,277 @@
+package roaring
+
+import (
+	"math/rand"
+	"testing"
+
+	"github.com/stretchr/testify/require"
+)
+
+// kwayKey picks the container key for the k-th key of an input: two inputs
+// built with the same k may land on different keys, and some keys are left
+// out, so key sets interleave and go missing the way real inputs do.
+func kwayKey(rng *rand.Rand, k int) (uint64, bool) {
+	return uint64(2*k+rng.Intn(2)) << 16, rng.Intn(4) != 0
+}
+
+// kwayBitmap builds a bitmap over up to keys keys holding random values at
+// the given density: arrays when sparse, bitmaps when dense.
+func kwayBitmap(rng *rand.Rand, keys int, density float64) *Bitmap {
+	bm := New()
+	for k := 0; k < keys; k++ {
+		base, present := kwayKey(rng, k)
+		if !present {
+			continue
+		}
+		for i := 0; i < int(density*65536); i++ {
+			bm.Add(uint32(base) + uint32(rng.Intn(65536)))
+		}
+	}
+	return bm
+}
+
+// kwayRuns builds a run-optimized bitmap of random ranges over up to keys keys.
+func kwayRuns(rng *rand.Rand, keys int) *Bitmap {
+	bm := New()
+	for k := 0; k < keys; k++ {
+		base, present := kwayKey(rng, k)
+		if !present {
+			continue
+		}
+		for start := 0; start < 65536; start += 1000 + rng.Intn(3000) {
+			bm.AddRange(base+uint64(start), base+uint64(min(start+1+rng.Intn(800), 65536)))
+		}
+	}
+	bm.RunOptimize()
+	return bm
+}
+
+// intervalBitmap holds one run container at key 0 made of the given
+// half-open intervals.
+func intervalBitmap(t *testing.T, intervals [][2]uint64) *Bitmap {
+	t.Helper()
+	bm := New()
+	for _, iv := range intervals {
+		bm.AddRange(iv[0], iv[1])
+	}
+	bm.RunOptimize()
+	require.Equal(t, uint64(1), bm.Stats().RunContainers, "%v", intervals)
+	return bm
+}
+
+func randomIntervals(rng *rand.Rand) [][2]uint64 {
+	var out [][2]uint64
+	pos := uint64(rng.Intn(64))
+	for n := 1 + rng.Intn(40); n > 0 && pos < 65536; n-- {
+		length := min(uint64(4+rng.Intn(3000)), 65536-pos)
+		out = append(out, [2]uint64{pos, pos + length})
+		pos += length + uint64(rng.Intn(200))
+	}
+	return out
+}
+
+// requireFastAnd checks FastAnd against the pairwise definition, that the
+// result is valid, and that it shares no storage with its inputs.
+func requireFastAnd(t *testing.T, inputs ...*Bitmap) {
+	t.Helper()
+	before := make([]*Bitmap, len(inputs))
+	for i, in := range inputs {
+		before[i] = in.Clone()
+	}
+	want := And(inputs[0], inputs[1])
+	for _, bm := range inputs[2:] {
+		want.And(bm)
+	}
+	got := FastAnd(inputs...)
+	require.NoError(t, got.Validate())
+	require.True(t, got.Equals(want), "%d values, want %d", got.GetCardinality(), want.GetCardinality())
+	got.Flip(0, 1<<21) // rewrites every container the inputs could share
+	for i, in := range inputs {
+		require.True(t, in.Equals(before[i]), "input %d was modified through the result", i)
+	}
+}
+
+func TestFastAndMatchesPairwise(t *testing.T) {
+	rng := rand.New(rand.NewSource(3))
+	densities := []float64{0.0005, 0.01, 0.05, 0.3, 0.8, 0.98}
+	for trial := 0; trial < 400; trial++ {
+		inputs := make([]*Bitmap, 3+rng.Intn(4))
+		for i := range inputs {
+			switch rng.Intn(5) {
+			case 0:
+				inputs[i] = kwayRuns(rng, 1+rng.Intn(4))
+			case 1:
+				inputs[i] = intervalBitmap(t, randomIntervals(rng))
+			default:
+				inputs[i] = kwayBitmap(rng, 1+rng.Intn(4), densities[rng.Intn(len(densities))])
+			}
+		}
+		requireFastAnd(t, inputs...)
+	}
+}
+
+func TestFastAndEdges(t *testing.T) {
+	rng := rand.New(rand.NewSource(12))
+	dense := func() *Bitmap { return kwayBitmap(rng, 1, 0.5) }
+	full := intervalBitmap(t, [][2]uint64{{0, 65536}})
+	var manyShort [][2]uint64
+	for p := uint64(0); p < 65536; p += 64 {
+		manyShort = append(manyShort, [2]uint64{p, p + 8})
+	}
+	for _, ivs := range [][][2]uint64{
+		{{0, 65536}},
+		{{777, 800}},
+		{{64, 128}, {1024, 4096}},
+		{{65472, 65536}},
+		{{60, 70}, {130, 140}, {65530, 65536}},
+		{{100, 200}, {201, 300}},
+		manyShort,
+	} {
+		mask := intervalBitmap(t, ivs)
+		requireFastAnd(t, dense(), mask, dense())             // one run masks the bitmaps
+		requireFastAnd(t, dense(), mask, full, dense())       // two runs fold into one mask
+		requireFastAnd(t, dense(), mask, full)                // one bitmap under two runs
+		requireFastAnd(t, mask, full, intervalBitmap(t, ivs)) // runs only
+	}
+	// Two runs whose fold, narrower than the first run, misses the one bitmap.
+	spot := New()
+	spot.Add(10)
+	spot.AddRange(5000, 20000)
+	requireFastAnd(t, intervalBitmap(t, [][2]uint64{{0, 100}}), intervalBitmap(t, [][2]uint64{{50, 100}}), spot)
+	// Ranges that each meet the first one but have no common value end the
+	// key before a bitmap is read; a range narrower than the mask's span
+	// clips a many-interval run.
+	requireFastAnd(t, intervalBitmap(t, [][2]uint64{{0, 40000}}), intervalBitmap(t, [][2]uint64{{30000, 65536}}), intervalBitmap(t, [][2]uint64{{0, 20000}}), dense())
+	requireFastAnd(t, intervalBitmap(t, [][2]uint64{{100, 40000}}), intervalBitmap(t, manyShort), dense(), dense())
+	requireFastAnd(t, intervalBitmap(t, [][2]uint64{{100, 40000}}), intervalBitmap(t, [][2]uint64{{50, 60}, {30000, 50000}}), dense())
+	// Two ranges that each meet the bitmap but overlap only where it is
+	// empty, on word boundaries: the one-bitmap span counts to zero.
+	gap := New()
+	gap.Add(0)
+	gap.AddRange(128, 5300)
+	requireFastAnd(t, gap, intervalBitmap(t, [][2]uint64{{0, 128}}), intervalBitmap(t, [][2]uint64{{64, 192}}))
+	// Disjoint runs end the call before a bitmap is read.
+	requireFastAnd(t, intervalBitmap(t, [][2]uint64{{0, 1000}}), intervalBitmap(t, [][2]uint64{{2000, 3000}}), dense())
+	// Keys present in some inputs only, and none shared at all.
+	x := kwayBitmap(rng, 6, 0.3)
+	requireFastAnd(t, x, AddOffset64(x, 1<<16), x)
+	requireFastAnd(t, x, AddOffset64(x, 40<<16), dense())
+	y := New()
+	y.AddMany([]uint32{1, 5, 9, 70000})
+	require.True(t, FastAnd(y, full, y).Equals(And(y, full)))
+	require.True(t, FastAnd(y, New(), full).IsEmpty())
+	require.True(t, FastAnd(full, full, full).Equals(full))
+	require.True(t, FastAnd(y, y, y).Equals(y))
+}
+
+// kwayTriple returns three dense inputs, every pair overlapping, all three
+// disjoint, so only the counted k-way pass can prove the result empty.
+func kwayTriple() (a, b, c *Bitmap) {
+	a, b, c = New(), New(), New()
+	for id := uint32(0); id < 4<<16; id++ {
+		switch id % 3 {
+		case 0:
+			a.Add(id)
+			b.Add(id)
+		case 1:
+			b.Add(id)
+			c.Add(id)
+		case 2:
+			a.Add(id)
+			c.Add(id)
+		}
+	}
+	return a, b, c
+}
+
+func TestFastAndEmptyAllocatesNoContainer(t *testing.T) {
+	a, b, c := kwayTriple()
+	require.True(t, a.Intersects(b) && b.Intersects(c) && a.Intersects(c))
+	require.True(t, FastAnd(a, b, c).IsEmpty())
+	// Only the answer bitmap; never a container.
+	require.LessOrEqual(t, testing.AllocsPerRun(50, func() { FastAnd(a, b, c) }), 1.0)
+
+	mask := intervalBitmap(t, [][2]uint64{{0, 65536}})
+	evens, odds := New(), New()
+	for v := uint32(0); v < 65536; v += 2 {
+		evens.Add(v)
+		odds.Add(v + 1)
+	}
+	require.True(t, FastAnd(mask, evens, odds).IsEmpty())
+	require.LessOrEqual(t, testing.AllocsPerRun(50, func() { FastAnd(mask, evens, odds) }), 1.0)
+}
+
+// shapeKey builds one container of the given kind and cardinality on key:
+// random values for arrays and bitmaps, card/8 short ranges for runs.
+func shapeKey(rng *rand.Rand, kind string, card int, key uint32) *Bitmap {
+	bm := New()
+	base := uint64(key) << 16
+	switch kind {
+	case "array", "bitmap":
+		for bm.GetCardinality() < uint64(card) {
+			bm.Add(uint32(base) + uint32(rng.Intn(65536)))
+		}
+	case "runs":
+		for i := 0; i < card/8; i++ {
+			bm.AddRange(base+uint64(i*16), base+uint64(i*16+8))
+		}
+		bm.RunOptimize()
+	}
+	return bm
+}
+
+func pairwiseAnd(bms []*Bitmap) *Bitmap {
+	r := And(bms[0], bms[1])
+	for _, bm := range bms[2:] {
+		r.And(bm)
+	}
+	return r
+}
+
+// BenchmarkFastAndShapes runs each kernel path once against the pairwise
+// chain on the same inputs: the smallest input an array, tiny and large,
+// against a bitmap, a run and a many-interval run; bitmaps under a run
+// mask; runs only; dense inputs that are pairwise non-empty and jointly
+// empty; and a term spanning eight keys.
+func BenchmarkFastAndShapes(b *testing.B) {
+	rng := rand.New(rand.NewSource(9))
+	window := New()
+	window.AddRange(0, 65536)
+	wide := New()
+	wide.AddRange(0, 8<<16)
+	spread := func(kind string, card int) *Bitmap {
+		bm := New()
+		for k := uint32(0); k < 8; k++ {
+			bm.Or(shapeKey(rng, kind, card, k))
+		}
+		return bm
+	}
+	x, y, z := kwayTriple()
+	for _, tc := range []struct {
+		name   string
+		inputs []*Bitmap
+	}{
+		{"array30-bitmap", []*Bitmap{window, shapeKey(rng, "array", 30, 0), shapeKey(rng, "bitmap", 33000, 0)}},
+		{"array4000-bitmap", []*Bitmap{window, shapeKey(rng, "array", 4000, 0), shapeKey(rng, "bitmap", 33000, 0)}},
+		{"array30-runs", []*Bitmap{window, shapeKey(rng, "array", 30, 0), shapeKey(rng, "runs", 4000, 0)}},
+		{"array4000-runs", []*Bitmap{window, shapeKey(rng, "array", 4000, 0), shapeKey(rng, "runs", 4000, 0)}},
+		{"bitmap-bitmap", []*Bitmap{window, shapeKey(rng, "bitmap", 33000, 0), shapeKey(rng, "bitmap", 33000, 0)}},
+		{"runs-bitmap", []*Bitmap{window, shapeKey(rng, "runs", 4000, 0), shapeKey(rng, "bitmap", 33000, 0)}},
+		{"runs-runs", []*Bitmap{window, shapeKey(rng, "runs", 4000, 0), shapeKey(rng, "runs", 30000, 0)}},
+		{"jointly-empty", []*Bitmap{x, y, z}},
+		{"eight-keys", []*Bitmap{wide, spread("array", 30), spread("bitmap", 33000)}},
+	} {
+		b.Run(tc.name+"/fastand", func(b *testing.B) {
+			b.ReportAllocs()
+			for b.Loop() {
+				_ = FastAnd(tc.inputs...)
+			}
+		})
+		b.Run(tc.name+"/pairwise", func(b *testing.B) {
+			b.ReportAllocs()
+			for b.Loop() {
+				_ = pairwiseAnd(tc.inputs)
+			}
+		})
+	}
+}
diff --git a/real_data_benchmark_test.go b/real_data_benchmark_test.go
index 7e47a323..f40a12a4 100644
--- a/real_data_benchmark_test.go
+++ b/real_data_benchmark_test.go
@@ -186,3 +186,9 @@ func BenchmarkRealDataFastOr(b *testing.B) {
 		return FastOr(bitmaps...).GetCardinality()
 	})
 }
+
+func BenchmarkRealDataFastAnd(b *testing.B) {
+	benchmarkRealDataAggregate(b, func(bitmaps []*Bitmap) uint64 {
+		return FastAnd(bitmaps...).GetCardinality()
+	})
+}
diff --git a/runbitmap_and_test.go b/runbitmap_and_test.go
new file mode 100644
index 00000000..33ffd333
--- /dev/null
+++ b/runbitmap_and_test.go
@@ -0,0 +1,260 @@
+package roaring
+
+import (
+	"math/rand"
+	"testing"
+
+	"github.com/stretchr/testify/require"
+)
+
+// randomRun builds a canonical run container: up to n sorted intervals of
+// at most maxLen values separated by gaps of at most maxGap, optionally
+// touching 0 and 65535.
+func randomRun(rng *rand.Rand, n, maxLen, maxGap int, touchEdges bool) *runContainer16 {
+	rc := &runContainer16{}
+	pos := 0
+	if !touchEdges {
+		pos = rng.Intn(100)
+	}
+	for i := 0; i < n && pos < 65536; i++ {
+		length := min(1+rng.Intn(1+rng.Intn(maxLen)), 65536-pos)
+		rc.iv = append(rc.iv, interval16{start: uint16(pos), length: uint16(length - 1)})
+		pos += length + 1 + rng.Intn(maxGap)
+	}
+	if last := int(rc.maximum()); touchEdges && last != 65535 {
+		if start := 65535 - rng.Intn(500); last+2 <= start {
+			rc.iv = append(rc.iv, interval16{start: uint16(start), length: uint16(65535 - start)})
+		}
+	}
+	return rc
+}
+
+// randomBitmapContainer ANDs ands random words per position, so the density
+// is 1/2^ands.
+func randomBitmapContainer(rng *rand.Rand, ands int) *bitmapContainer {
+	bc := newBitmapContainer()
+	for i := range bc.bitmap {
+		w := rng.Uint64()
+		for k := 1; k < ands; k++ {
+			w &= rng.Uint64()
+		}
+		bc.bitmap[i] = w
+	}
+	bc.computeCardinality()
+	return bc
+}
+
+func requireAndResult(t *testing.T, want, got container) {
+	t.Helper()
+	require.True(t, want.equals(got), "want %d values, got %d", want.getCardinality(), got.getCardinality())
+	_, isArray := got.(*arrayContainer)
+	require.Equal(t, got.getCardinality() <= arrayDefaultMaxSize, isArray, "container kind must follow its cardinality")
+}
+
+// The reference is the old path: materialize the run, then the unchanged
+// bitmap-bitmap intersection.
+func TestRunAndBitmapContainer(t *testing.T) {
+	rng := rand.New(rand.NewSource(1))
+	for trial := 0; trial < 300; trial++ {
+		rc := randomRun(rng, 1+rng.Intn(40), 3000, 2000, trial%5 == 0)
+		bc := randomBitmapContainer(rng, 1+trial%4)
+		want := newBitmapContainerFromRun(rc).andBitmap(bc)
+		requireAndResult(t, want, rc.and(bc))
+		requireAndResult(t, want, bc.and(rc))
+		requireAndResult(t, want, bc.clone().iand(rc))
+	}
+}
+
+// Runs with more than runAndScratchIntervals intervals take the scratch
+// path; check it on regular and random runs against every density.
+func TestRunAndBitmapContainerManyIntervals(t *testing.T) {
+	rng := rand.New(rand.NewSource(21))
+	regular := func(stride, length int) *runContainer16 {
+		rc := &runContainer16{}
+		for start := 0; start+length <= 65536; start += stride {
+			rc.iv = append(rc.iv, interval16{start: uint16(start), length: uint16(length - 1)})
+		}
+		return rc
+	}
+	runs := []*runContainer16{regular(32, 8), regular(64, 1), regular(3, 2), regular(200, 100)}
+	for i := 0; i < 60; i++ {
+		runs = append(runs, randomRun(rng, 65+rng.Intn(400), 60, 120, false))
+	}
+	for i, rc := range runs {
+		require.Greater(t, len(rc.iv), runAndScratchIntervals, "run %d", i)
+		for _, bc := range []*bitmapContainer{randomBitmapContainer(rng, 6), randomBitmapContainer(rng, 2), randomBitmapContainer(rng, 1), newBitmapContainerwithRange(0, 65535)} {
+			want := newBitmapContainerFromRun(rc).andBitmap(bc)
+			requireAndResult(t, want, rc.and(bc))
+			requireAndResult(t, want, bc.and(rc))
+			requireAndResult(t, want, bc.clone().iand(rc))
+			require.Equal(t, !want.isEmpty(), rc.intersects(bc), "run %d", i)
+		}
+	}
+}
+
+func TestRunAndBitmapContainerEdges(t *testing.T) {
+	empty := newBitmapContainer()
+	full := newBitmapContainerwithRange(0, 65535)
+	single := &runContainer16{iv: []interval16{{start: 100, length: 9}}}
+	edges := &runContainer16{iv: []interval16{{start: 0, length: 0}, {start: 65535, length: 0}}}
+	atThreshold := &runContainer16{iv: []interval16{{start: 0, length: arrayDefaultMaxSize - 1}}}
+	pastThreshold := &runContainer16{iv: []interval16{{start: 0, length: arrayDefaultMaxSize}}}
+	for _, rc := range []*runContainer16{single, edges, atThreshold, pastThreshold} {
+		for _, bc := range []*bitmapContainer{empty, full} {
+			want := newBitmapContainerFromRun(rc).andBitmap(bc)
+			requireAndResult(t, want, rc.and(bc))
+			requireAndResult(t, want, bc.clone().iand(rc))
+		}
+	}
+}
+
+// An empty intersection allocates only its empty array container; the old
+// path allocated an 8 KiB bitmap to hold the run first.
+func TestRunAndBitmapContainerEmptyAllocation(t *testing.T) {
+	rc := &runContainer16{iv: []interval16{{start: 0, length: 999}, {start: 30000, length: 999}}}
+	bc := newBitmapContainerwithRange(2000, 20000)
+	require.LessOrEqual(t, testing.AllocsPerRun(100, func() {
+		if !rc.and(bc).isEmpty() {
+			t.Fatal("expected empty")
+		}
+	}), 1.0)
+	// The clone costs two allocations, the empty result one.
+	require.LessOrEqual(t, testing.AllocsPerRun(100, func() {
+		_ = bc.clone().iand(rc)
+	}), 3.0)
+}
+
+// The reference is the old definition, !rc.and(c).isEmpty().
+func TestRunContainerIntersects(t *testing.T) {
+	rng := rand.New(rand.NewSource(11))
+	for trial := 0; trial < 400; trial++ {
+		rc := randomRun(rng, 1+rng.Intn(30), 3000, 2000, trial%7 == 0)
+		ac := newArrayContainer()
+		for i := 0; i < rng.Intn(300); i++ {
+			ac.iadd(uint16(rng.Intn(65536)))
+		}
+		others := []container{randomBitmapContainer(rng, 1+trial%3*3), ac, randomRun(rng, 1+rng.Intn(30), 3000, 2000, trial%5 == 0)}
+		for _, c := range others {
+			want := !rc.and(c).isEmpty()
+			require.Equal(t, want, rc.intersects(c), "trial %d run vs %T", trial, c)
+			require.Equal(t, want, c.intersects(rc), "trial %d %T vs run", trial, c)
+		}
+	}
+	// A few values against many intervals take the search branch; the
+	// reference is the same definition.
+	many := randomRun(rng, 300, 60, 120, false)
+	require.Greater(t, len(many.iv), 200)
+	for trial := 0; trial < 200; trial++ {
+		few := newArrayContainer()
+		for i := 0; i < 1+rng.Intn(6); i++ {
+			few.iadd(uint16(rng.Intn(65536)))
+		}
+		want := !many.and(few).isEmpty()
+		require.Equal(t, want, many.intersects(few), "trial %d few vs many", trial)
+		require.Equal(t, want, few.intersects(many), "trial %d many vs few", trial)
+	}
+	// Edges: empty array, touching intervals, single values at both ends.
+	empty := newArrayContainer()
+	edges := &runContainer16{iv: []interval16{{start: 0, length: 0}, {start: 65535, length: 0}}}
+	require.False(t, edges.intersects(empty))
+	one := newArrayContainer()
+	one.iadd(65535)
+	require.True(t, edges.intersects(one))
+	touching := &runContainer16{iv: []interval16{{start: 1, length: 65533}}}
+	require.False(t, edges.intersects(touching))
+	require.True(t, edges.intersects(newBitmapContainerwithRange(0, 0)))
+}
+
+func TestRunContainerIntersectsAllocatesNothing(t *testing.T) {
+	rc := &runContainer16{iv: []interval16{{start: 0, length: 999}, {start: 30000, length: 999}}}
+	bc := newBitmapContainerwithRange(2000, 20000)
+	ac := newArrayContainerRange(2000, 3000)
+	other := &runContainer16{iv: []interval16{{start: 5000, length: 100}}}
+	for _, c := range []container{bc, ac, other} {
+		require.Zero(t, testing.AllocsPerRun(100, func() { _ = rc.intersects(c) }), "%T", c)
+	}
+}
+
+// runBitmap holds one run container at key 0 covering n intervals of the
+// given length, spaced stride apart.
+func runBitmap(n, length, stride int) *Bitmap {
+	bm := New()
+	for i := 0; i < n; i++ {
+		bm.AddRange(uint64(i*stride), uint64(i*stride+length))
+	}
+	bm.RunOptimize()
+	return bm
+}
+
+// BenchmarkRunAndBitmap intersects a run-optimized bitmap with a dense one
+// over a single key through the public API, one case per result kind, plus
+// the in-place form and the Intersects test.
+func BenchmarkRunAndBitmap(b *testing.B) {
+	rng := rand.New(rand.NewSource(7))
+	short := runBitmap(20, 300, 3000)  // 6000 values
+	long := runBitmap(4, 10000, 16000) // 40000 values
+	many := runBitmap(1024, 8, 64)     // 8192 values in 1024 intervals
+	gaps, quarter, half, dense := New(), New(), New(), New()
+	for v := uint32(0); v < 65536; v++ {
+		if !short.Contains(v) && rng.Intn(2) == 0 {
+			gaps.Add(v)
+		}
+		if rng.Intn(4) == 0 {
+			quarter.Add(v)
+		}
+		if rng.Intn(2) == 0 {
+			half.Add(v)
+		}
+		if rng.Intn(10) != 0 {
+			dense.Add(v)
+		}
+	}
+	for _, tc := range []struct {
+		name string
+		x, y *Bitmap
+	}{
+		{"and/empty", short, gaps},
+		{"and/array", short, half},
+		{"and/bitmap", long, dense},
+		{"and/many-intervals/array", many, quarter},
+		{"and/many-intervals/bitmap", many, dense},
+	} {
+		b.Run(tc.name, func(b *testing.B) {
+			b.ReportAllocs()
+			for b.Loop() {
+				_ = And(tc.x, tc.y)
+			}
+		})
+	}
+	for _, tc := range []struct {
+		name string
+		x    *Bitmap
+	}{
+		{"iand/bitmap", long},
+		{"iand/many-intervals", many},
+	} {
+		b.Run(tc.name, func(b *testing.B) {
+			b.ReportAllocs()
+			x := dense.Clone()
+			for b.Loop() {
+				x.And(tc.x)
+			}
+		})
+	}
+	few := New()
+	few.AddMany([]uint32{100, 30000, 60000})
+	for _, tc := range []struct {
+		name string
+		x, y *Bitmap
+	}{
+		{"intersects/disjoint", short, gaps},
+		{"intersects/overlapping", short, half},
+		{"intersects/few-values-many-intervals", many, few},
+	} {
+		b.Run(tc.name, func(b *testing.B) {
+			for b.Loop() {
+				_ = tc.x.Intersects(tc.y)
+			}
+		})
+	}
+}
diff --git a/runcontainer.go b/runcontainer.go
index 7c369b75..5cf35f24 100644
--- a/runcontainer.go
+++ b/runcontainer.go
@@ -1855,10 +1855,42 @@ func (rc *runContainer16) andCardinality(a container) int {
 	panic("unsupported container type")
 }
 
-// andBitmapContainer finds the intersection of rc and b.
+// Above this many intervals, building the intersection in scratch beats
+// counting it interval by interval first (measured).
+const runAndScratchIntervals = 64
+
+// andBitmapContainer intersects rc with bc, counting before allocating.
 func (rc *runContainer16) andBitmapContainer(bc *bitmapContainer) container {
-	bc2 := newBitmapContainerFromRun(rc)
-	return bc2.andBitmap(bc)
+	if len(rc.iv) > runAndScratchIntervals {
+		return rc.andBitmapContainerScratch(bc)
+	}
+	card := rc.andBitmapContainerCardinality(bc)
+	if card > arrayDefaultMaxSize {
+		answer := newBitmapContainer()
+		for i := range rc.iv {
+			orBitmapRange(answer.bitmap, bc.bitmap, int(rc.iv[i].start), int(rc.iv[i].last())+1)
+		}
+		answer.cardinality = card
+		return answer
+	}
+	answer := newArrayContainerCapacity(card)
+	for i := range rc.iv {
+		answer.content = appendBitmapRange(answer.content, bc.bitmap, int(rc.iv[i].start), int(rc.iv[i].last())+1)
+	}
+	return answer
+}
+
+func (rc *runContainer16) andBitmapContainerScratch(bc *bitmapContainer) container {
+	first, last := int(rc.iv[0].start)/64, int(rc.maximum())/64
+	var scratch [bitmapContainerSize]uint64
+	for i := range rc.iv {
+		orBitmapRange(scratch[:], bc.bitmap, int(rc.iv[i].start), int(rc.iv[i].last())+1)
+	}
+	card := int(popcntSlice(scratch[first : last+1]))
+	if card == 0 {
+		return newArrayContainerCapacity(0)
+	}
+	return containerFromWords(scratch[:], first, last, card)
 }
 
 func (rc *runContainer16) andArrayCardinality(ac *arrayContainer) int {
@@ -2418,9 +2450,53 @@ func (rc *runContainer16) lazyOR(a container) container {
 }
 
 func (rc *runContainer16) intersects(a container) bool {
-	// TODO: optimize by doing inplace/less allocation
-	isect := rc.and(a)
-	return !isect.isEmpty()
+	switch c := a.(type) {
+	case *bitmapContainer:
+		for i := range rc.iv {
+			if c.intersectsRange(uint(rc.iv[i].start), uint(rc.iv[i].last())+1) {
+				return true
+			}
+		}
+		return false
+	case *arrayContainer:
+		if len(c.content)*bits.Len(uint(len(rc.iv))) < len(rc.iv) {
+			// Few values against many intervals: a search each beats
+			// walking every interval.
+			for _, v := range c.content {
+				if rc.contains(v) {
+					return true
+				}
+			}
+			return false
+		}
+		pos := 0
+		for i := range rc.iv {
+			// advanceUntil searches from pos+1; the value at pos may still match.
+			pos = advanceUntil(c.content, pos-1, len(c.content), rc.iv[i].start)
+			if pos == len(c.content) {
+				return false
+			}
+			if c.content[pos] <= rc.iv[i].last() {
+				return true
+			}
+		}
+		return false
+	case *runContainer16:
+		i, j := 0, 0
+		for i < len(rc.iv) && j < len(c.iv) {
+			x, y := rc.iv[i], c.iv[j]
+			if haveOverlap16(x, y) {
+				return true
+			}
+			if x.last() < y.last() {
+				i++
+			} else {
+				j++
+			}
+		}
+		return false
+	}
+	panic("unsupported container type")
 }
 
 func (rc *runContainer16) xor(a container) container {
diff --git a/util.go b/util.go
index e727d6e0..6ed3934e 100644
--- a/util.go
+++ b/util.go
@@ -143,6 +143,70 @@ func flipBitmapRange(bitmap []uint64, start int, end int) {
 	bitmap[endword] ^= ^uint64(0) >> (uint(-end) % 64)
 }
 
+// orBitmapRange ORs the bits of src within [start, end) into dst.
+func orBitmapRange(dst, src []uint64, start int, end int) {
+	if start >= end {
+		return
+	}
+	firstword := start / 64
+	endword := (end - 1) / 64
+	if firstword == endword {
+		dst[firstword] |= src[firstword] & (^uint64(0) << uint(start%64)) & (^uint64(0) >> (uint(-end) % 64))
+		return
+	}
+	dst[firstword] |= src[firstword] & (^uint64(0) << uint(start%64))
+	for i := firstword + 1; i < endword; i++ {
+		dst[i] |= src[i]
+	}
+	dst[endword] |= src[endword] & (^uint64(0) >> (uint(-end) % 64))
+}
+
+// appendBitmapRange appends the positions of the bits set in [start, end).
+func appendBitmapRange(out []uint16, bitmap []uint64, start int, end int) []uint16 {
+	if start >= end {
+		return out
+	}
+	firstword := start / 64
+	endword := (end - 1) / 64
+	for i := firstword; i <= endword; i++ {
+		w := bitmap[i]
+		if i == firstword {
+			w &= ^uint64(0) << uint(start%64)
+		}
+		if i == endword {
+			w &= ^uint64(0) >> (uint(-end) % 64)
+		}
+		for w != 0 {
+			out = append(out, uint16(i*64+bits.TrailingZeros64(w)))
+			w &= w - 1
+		}
+	}
+	return out
+}
+
+// clearBitmapGaps clears the bits between from, the sorted intervals iv, and to.
+func clearBitmapGaps(bitmap []uint64, iv []interval16, from, to int) {
+	for i := range iv {
+		resetBitmapRange(bitmap, from, int(iv[i].start))
+		from = int(iv[i].last()) + 1
+	}
+	resetBitmapRange(bitmap, from, to)
+}
+
+// containerFromWords builds an array or bitmap container from the card bits
+// set in words [first, last] of w.
+func containerFromWords(w []uint64, first, last, card int) container {
+	if card <= arrayDefaultMaxSize {
+		ac := newArrayContainerCapacity(card)
+		ac.content = appendBitmapRange(ac.content, w, first*64, (last+1)*64)
+		return ac
+	}
+	bc := newBitmapContainer()
+	copy(bc.bitmap[first:last+1], w[first:last+1])
+	bc.cardinality = card
+	return bc
+}
+
 func resetBitmapRange(bitmap []uint64, start int, end int) {
 	if start >= end {
 		return