diff --git a/roaring/roaring.go b/roaring/roaring.go
index ead9e6c43..4b9139508 100644
--- a/roaring/roaring.go
+++ b/roaring/roaring.go
@@ -2198,6 +2198,95 @@ func (c *Container) Contains(v uint16) bool {
 	}
 }
 
+// BitwiseCompare reports whether two containers are equal. It returns
+// an error describing any difference it finds. This is mostly intended
+// for use in tests that expect equality.
+func (c *Container) BitwiseCompare(c2 *Container) error {
+	if c.N() != c2.N() {
+		return errors.New("containers are different lengths")
+	}
+	if c.N() == 0 {
+		return nil
+	}
+	switch typePair(c.typ(), c2.typ()) {
+	case typePair(containerArray, containerArray):
+		return compareArrayArray(c.array(), c2.array())
+	case typePair(containerArray, containerBitmap):
+		return compareArrayBitmap(c.array(), c2.bitmap())
+	case typePair(containerBitmap, containerArray):
+		return compareArrayBitmap(c2.array(), c.bitmap())
+	case typePair(containerArray, containerRun):
+		return compareArrayRuns(c.array(), c2.runs())
+	case typePair(containerRun, containerArray):
+		return compareArrayRuns(c2.array(), c.runs())
+	default:
+		c3 := xor(c, c2)
+		if c3.N() != 0 {
+			return fmt.Errorf("%d bits differenct between containers", c3.N())
+		}
+	}
+	return nil
+}
+
+func typePair(ct1, ct2 byte) int {
+	return int((ct1 << 4) | ct2)
+}
+
+func compareArrayArray(a1, a2 []uint16) error {
+	if len(a1) != len(a2) {
+		return fmt.Errorf("unexpected length mismatch, %d vs %d", len(a1), len(a2))
+	}
+	for i := range a1 {
+		if a1[i] != a2[i] {
+			return fmt.Errorf("item %d: %d vs %d", i, a1[i], a2[i])
+		}
+	}
+	return nil
+}
+
+// compareArrayRuns determines whether an array matches a provided
+// set of runs. As with compareArrayBitmap, it only verifies presence
+// of the array's values in the run collection. the run collection
+// can't be empty; if it were, N would have been 0, and we wouldn't
+// have gotten here.
+func compareArrayRuns(a []uint16, r []interval16) error {
+	ri := 0
+	ru := r[ri]
+	ri++
+	for _, v := range a {
+		if v < ru.start {
+			return fmt.Errorf("value %d missing", v)
+		}
+		if v > ru.last {
+			if ri >= len(r) {
+				return fmt.Errorf("value %d missing", v)
+			}
+			ru = r[ri]
+			ri++
+			// if they're identical, the array value must be
+			// the start of the next run.
+			if v != ru.start {
+				return fmt.Errorf("value %d missing", v)
+			}
+		}
+	}
+	return nil
+}
+
+// compareArrayBitmap actually only verifies that everything in the array
+// is in the bitmap. It's used only after comparing the N for the containers,
+// so if there's anything in the bitmap that's not in the array, either there's
+// something in the array that's not in the bitmap, or we didn't get here.
+func compareArrayBitmap(a []uint16, b []uint64) error {
+	for _, v := range a {
+		w, bit := b[v>>6], v&63
+		if w>>bit&1 == 0 {
+			return fmt.Errorf("value %d missing", v)
+		}
+	}
+	return nil
+}
+
 func (c *Container) bitmapCountRuns() (r int32) {
 	return bitmapCountRuns(c.bitmap())
 }
@@ -2342,8 +2431,7 @@ func (c *Container) unionInPlace(other *Container) *Container {
 			c = c.runToBitmap()
 			return unionBitmapArrayInPlace(c, other)
 		case containerRun:
-			c = c.runToBitmap()
-			return unionBitmapRunInPlace(c, other)
+			return unionRunRunInPlace(c, other)
 		}
 	}
 	if roaringParanoia {
@@ -2458,7 +2546,6 @@ func (c *Container) runRemove(v uint16) (*Container, bool) {
 		runs = append(runs, interval16{})
 		copy(runs[i+2:], runs[i+1:])
 		runs[i+1] = interval16{start: v + 1, last: last}
-		// runs = append(runs[:i+1], append([]interval16{{start: v + 1, last: last}}, runs[i+1:]...)...)
 	}
 	c.setN(c.N() - 1)
 	c.setRuns(runs)
@@ -3838,6 +3925,197 @@ func unionBitmapBitmapInPlace(a, b *Container) *Container {
 	return a
 }
 
+// unionRunRunInPlace unions run b into run a, mutating a in place.
+func unionRunRunInPlace(a, b *Container) *Container {
+	a = a.Thaw()
+	runs, n := unionInterval16InPlace(a.runs(), b.runs())
+
+	a.setRuns(runs)
+	a.setN(n)
+	return a
+}
+
+// unioninterval16InPlace merges two slices in place (in a).
+// The main concept is to go value by value (instead of interval by interval)
+// and count `.start` and `.last` points.
+// If we get the `state == 0` it means we just built a new interval (`val`),
+// and we can set it in `a` at the possition `off`
+func unionInterval16InPlace(a, b []interval16) ([]interval16, int32) {
+	n := int32(0)
+	an, bn := len(a), len(b)
+
+	var (
+		// ai - index of a, aii - subindex (0: a[ai].start, 1: a[ai].last).
+		ai, aii int = 0, 0
+
+		// bi - index of b, bii - subindex (0: b[bi].start, 1: b[bi].last).
+		bi, bii int = 0, 0
+
+		// Offset of a - next available index to set.
+		off int = 0
+		// Value to set/append to a at off
+		val interval16
+
+		// Current state - state equals 0 means we are clear (out of intervals)
+		// When we start a new interval we add +1,
+		// when we get out of interval we add  -1.
+		state int
+
+		// mapping: subindex (ii) to state
+		// .start: [0] ->  1
+		// .last:  [1] -> -1
+		iiMap = [2]int{1, -1}
+
+		// If fromB is equal 2 it means that both val.start and val.last come from b,
+		// so we need to extend a, first
+		fromB int8
+
+		// eval functions evaluates global state and value
+		eval = func(arr [2]uint16, ii int, onlyB bool) {
+			if state == 0 && ii == 0 {
+				// we are clear and start a new interval
+				val.start = arr[ii]
+				if onlyB {
+					fromB++
+				}
+			}
+
+			state += iiMap[ii]
+
+			if state == 0 {
+				// we just got out of interval
+				// ii == 1
+				val.last = arr[ii]
+				if onlyB {
+					fromB++
+				}
+			}
+		}
+		// eval2 function is a special variant for eval function
+		// it's only used when two interval endings are equal, e.g.:
+		// a: ------------------|
+		// b:        -----------|
+		// the most important part is to change the global for both endings
+		// before we check if we're getting out of interval and start the new one.
+		eval2 = func(arr [2]uint16, i1, i2 int) {
+			if state == 0 && (i1 == 0 || i2 == 0) {
+				// we are clear and start a new interval
+				val.start = arr[i1]
+
+			}
+
+			state += iiMap[i1]
+			state += iiMap[i2]
+
+			if state == 0 {
+				// (i1 == 1 || i2 == 1)
+				// we just got out of interval
+				val.last = arr[i1]
+			}
+		}
+	)
+
+	for {
+		// av, bv reflects a[ai] and b[bi] intervals as an array,
+		// so we can internally iterate over values (points).
+		var av, bv [2]uint16
+
+		if ai < an && bi < bn {
+			av[0], av[1] = a[ai].start, a[ai].last
+			bv[0], bv[1] = b[bi].start, b[bi].last
+
+			if av[aii] < bv[bii] {
+				// a: |-------------------
+				// b:           |-------------------
+
+				eval(av, aii, false)
+				aii++
+			} else if av[aii] == bv[bii] {
+				// a: |-------------------
+				// b: |-------------------
+				// or
+				// a: ------------------|
+				// b:                   |-------------------
+				// or
+				// a: ------------------|
+				// b:      |------------|
+				// ...
+
+				eval2(av, aii, bii)
+				aii++
+				bii++
+			} else { // bv[bii] < av[aii]
+				// a:           |-------------------
+				// b: |-------------------
+
+				eval(bv, bii, true)
+				bii++
+			}
+		} else if ai < an { // only a left
+			av[0], av[1] = a[ai].start, a[ai].last
+			eval(av, aii, false)
+			aii++
+		} else if bi < bn { // only b left
+			bv[0], bv[1] = b[bi].start, b[bi].last
+			eval(bv, bii, false)
+			bii++
+		} else {
+			break
+		}
+
+		if state == 0 {
+			if fromB == 2 {
+				// val.start and val.last come from b, so we need to extend a, first
+				a = append(a, interval16{})
+				copy(a[off+1:], a[off:])
+				ai++
+				an++
+			}
+			fromB = 0
+			a, off = appendinterval16At(a, val, off)
+			n += int32(val.last) - int32(val.start) + 1
+		}
+
+		if aii == 2 {
+			// move to the next a's interval
+			aii = 0
+			ai++
+		}
+
+		if bii == 2 {
+			// move to the next b's interval
+			bii = 0
+			bi++
+		}
+	}
+
+	if len(a) > 0 {
+		a = a[:off]
+	}
+	return a, n
+}
+
+// appendinterval16At appends or sets val in a at off position
+// The function returns modified a ([]interval16) and new offset (off)
+func appendinterval16At(a []interval16, val interval16, off int) ([]interval16, int) {
+
+	if off > 0 && int32(val.start)-int32(a[off-1].last) <= 1 {
+		a[off-1].last = val.last
+		return a, off
+	}
+
+	if off == len(a) {
+		a = append(a, val)
+		off++
+		return a, off
+	}
+
+	a[off] = val
+	off++
+
+	return a, off
+}
+
 func difference(a, b *Container) *Container {
 	if a.N() == 0 || b.N() == maxContainerVal+1 {
 		return nil
diff --git a/roaring/roaring_internal_test.go b/roaring/roaring_internal_test.go
index 27fbe8f29..81d6515d8 100644
--- a/roaring/roaring_internal_test.go
+++ b/roaring/roaring_internal_test.go
@@ -883,6 +883,221 @@ func TestUnionRunRun(t *testing.T) {
 	}
 }
 
+func TestUnionInterval16InPlace(t *testing.T) {
+	tests := []struct {
+		name      string
+		a         []interval16
+		b         []interval16
+		expected  []interval16
+		expectedN int32
+	}{
+		{
+			name:      "firstBitUnset lastBitSet",
+			a:         []interval16{interval16{1, 10}},
+			b:         []interval16{interval16{10, 10}},
+			expected:  []interval16{interval16{1, 10}},
+			expectedN: 10,
+		},
+		{
+			name:      "single overlap",
+			a:         []interval16{interval16{1, 10}, interval16{21, 28}},
+			b:         []interval16{interval16{8, 12}},
+			expected:  []interval16{interval16{1, 12}, interval16{21, 28}},
+			expectedN: 20,
+		},
+		{
+			name:      "nested intervals",
+			a:         []interval16{interval16{3, 13}, interval16{17, 20}},
+			b:         []interval16{interval16{1, 4}, interval16{6, 7}, interval16{8, 9}, interval16{10, 11}, interval16{14, 17}},
+			expected:  []interval16{interval16{1, 20}},
+			expectedN: 20,
+		},
+		{
+			name:      "no overlap",
+			a:         []interval16{interval16{3, 4}, interval16{7, 8}},
+			b:         []interval16{interval16{1, 2}, interval16{5, 6}, interval16{9, 10}},
+			expected:  []interval16{interval16{1, 10}},
+			expectedN: 10,
+		},
+		{
+			name:      "b in a",
+			a:         []interval16{interval16{1, 10}},
+			b:         []interval16{interval16{5, 7}},
+			expected:  []interval16{interval16{1, 10}},
+			expectedN: 10,
+		},
+		{
+			name:      "a eq b",
+			a:         []interval16{interval16{1, 10}},
+			b:         []interval16{interval16{1, 10}},
+			expected:  []interval16{interval16{1, 10}},
+			expectedN: 10,
+		},
+		{
+			name:      "a in b",
+			a:         []interval16{interval16{5, 7}},
+			b:         []interval16{interval16{1, 10}},
+			expected:  []interval16{interval16{1, 10}},
+			expectedN: 10,
+		},
+		{
+			name:      "a ahead b",
+			a:         []interval16{interval16{1, 2}, interval16{3, 4}, interval16{5, 7}},
+			b:         []interval16{interval16{10, 11}, interval16{12, 13}, interval16{14, 15}},
+			expected:  []interval16{interval16{1, 7}, interval16{10, 15}},
+			expectedN: 13,
+		},
+		{
+			name:      "b ahead a",
+			a:         []interval16{interval16{10, 11}, interval16{12, 13}, interval16{14, 15}},
+			b:         []interval16{interval16{1, 2}, interval16{3, 4}, interval16{5, 7}},
+			expected:  []interval16{interval16{1, 7}, interval16{10, 15}},
+			expectedN: 13,
+		},
+		{
+			name:      "empty a and b",
+			a:         []interval16{},
+			b:         []interval16{},
+			expected:  []interval16{},
+			expectedN: 0,
+		},
+		{
+			name:      "empty a",
+			a:         []interval16{},
+			b:         []interval16{interval16{1, 2}, interval16{3, 4}, interval16{5, 7}},
+			expected:  []interval16{interval16{1, 7}},
+			expectedN: 7,
+		},
+		{
+			name:      "empty b",
+			a:         []interval16{interval16{1, 2}, interval16{3, 4}, interval16{5, 7}},
+			b:         []interval16{},
+			expected:  []interval16{interval16{1, 7}},
+			expectedN: 7,
+		},
+		{
+			name:      "single a",
+			a:         []interval16{interval16{1, 2}},
+			b:         []interval16{},
+			expected:  []interval16{interval16{1, 2}},
+			expectedN: 2,
+		},
+		{
+			name:      "single b",
+			a:         []interval16{},
+			b:         []interval16{interval16{1, 2}},
+			expected:  []interval16{interval16{1, 2}},
+			expectedN: 2,
+		},
+		{
+			name:      "single a single b",
+			a:         []interval16{interval16{3, 4}},
+			b:         []interval16{interval16{1, 2}},
+			expected:  []interval16{interval16{1, 4}},
+			expectedN: 4,
+		},
+		{
+			name:      "oddBitsSet lastBitUnset",
+			a:         []interval16{interval16{1, 1}, interval16{3, 3}, interval16{5, 5}},
+			b:         []interval16{interval16{0, 4}},
+			expected:  []interval16{interval16{0, 5}},
+			expectedN: 6,
+		},
+		{
+			name:      "all bits",
+			a:         []interval16{interval16{1, 1}, interval16{3, 3}, interval16{5, 5}},
+			b:         []interval16{interval16{0, 0}, interval16{2, 2}, interval16{4, 4}},
+			expected:  []interval16{interval16{0, 5}},
+			expectedN: 6,
+		},
+		{
+			name:      "short a long b",
+			a:         []interval16{interval16{5, 5}, interval16{7, 7}, interval16{9, 10}, interval16{12, 12}, interval16{15, 17}, interval16{19, 20}},
+			b:         []interval16{interval16{1, 10}, interval16{12, 12}, interval16{14, 18}},
+			expected:  []interval16{interval16{1, 10}, interval16{12, 12}, interval16{14, 20}},
+			expectedN: 18,
+		},
+		{
+			name:      "common endings",
+			a:         []interval16{interval16{1, 5}, interval16{15, 20}, interval16{25, 35}},
+			b:         []interval16{interval16{1, 10}, interval16{15, 20}, interval16{30, 35}},
+			expected:  []interval16{interval16{1, 10}, interval16{15, 20}, interval16{25, 35}},
+			expectedN: 27,
+		},
+		{
+			name:      "common endings and overlap",
+			a:         []interval16{interval16{1, 5}, interval16{10, 15}},
+			b:         []interval16{interval16{5, 10}, interval16{12, 17}},
+			expected:  []interval16{interval16{1, 17}},
+			expectedN: 17,
+		},
+		{
+			name:      "no common endings and overlap",
+			a:         []interval16{interval16{5, 10}, interval16{12, 17}},
+			b:         []interval16{interval16{0, 11}, interval16{15, 20}},
+			expected:  []interval16{interval16{0, 20}},
+			expectedN: 21,
+		},
+	}
+
+	for _, tc := range tests {
+		t.Run(tc.name, func(t *testing.T) {
+			bb := make([]interval16, len(tc.b))
+			copy(bb, tc.b)
+
+			runs, n := unionInterval16InPlace(tc.a, tc.b)
+
+			for i, v := range tc.expected {
+				if runs[i] != v {
+					t.Fatalf("runs expected: %+v, got: %+v", tc.expected, runs)
+				}
+			}
+			if n != tc.expectedN {
+				t.Fatalf("N expected: %d, got: %d", tc.expectedN, n)
+			}
+
+			for i, v := range bb {
+				if tc.b[i] != v {
+					t.Fatalf("b changed - runs expected: %+v, got: %+v", bb, tc.b)
+				}
+			}
+		})
+	}
+}
+
+func TestUnionRunRunInPlaceBitwiseCompare(t *testing.T) {
+	runs := []struct {
+		name string
+		run  []interval16
+	}{
+		{name: "FirstBitSet", run: runFirstBitSet()},
+		{name: "LastBitSet", run: runLastBitSet()},
+		{name: "FirstBitUnset", run: runFirstBitUnset()},
+		{name: "LastBitUnset", run: runLastBitUnset()},
+		{name: "InnerBitsSet", run: runInnerBitsSet()},
+		{name: "OuterBitsSet", run: runOuterBitsSet()},
+		{name: "OddBitsSet", run: runOddBitsSet()},
+		{name: "EvenBitsSet", run: runEvenBitsSet()},
+	}
+
+	for _, a := range runs {
+		for _, b := range runs {
+			t.Run(a.name+"-"+b.name, func(t *testing.T) {
+				arun := doContainer(containerRun, a.run)
+				brun := doContainer(containerRun, b.run)
+
+				out1 := unionBitmapRunInPlace(arun.runToBitmap(), brun)
+				out2 := unionRunRunInPlace(arun, brun)
+
+				err := out1.BitwiseCompare(out2.runToBitmap())
+				if err != nil {
+					t.Fatal(err)
+				}
+			})
+		}
+	}
+}
+
 func TestUnionArrayRun(t *testing.T) {
 	a := NewContainerArray(nil)
 	b := NewContainerRun(nil)
@@ -1441,11 +1656,13 @@ func TestDifferenceRunArray(t *testing.T) {
 		}
 	}
 }
+
 func MakeBitmap(start []uint64) []uint64 {
 	b := make([]uint64, bitmapN)
 	copy(b, start)
 	return b
 }
+
 func MakeLastBitSet() []uint64 {
 	obj := NewFileBitmap(65535)
 	c := obj.container(0)
@@ -2453,6 +2670,7 @@ func TestRunBinSearch(t *testing.T) {
 		}
 	}
 }
+
 func TestBitmap_RemoveEmptyContainers(t *testing.T) {
 	bm1 := NewFileBitmap(1<<16, 2<<16, 3<<16)
 	bm2 := NewFileBitmap(1<<16, 2<<16+1, 3<<16)
@@ -3488,28 +3706,6 @@ func newTestBitmapContainer() *Container {
 	return NewContainerBitmap(0, nil)
 }
 
-/*
-// This function exercises an arcane edge case in dead code.
-// It doesn't need to be run right now.
-func TestEquals(t *testing.T) {
-	bma := NewBitmap()
-	bmr := NewBitmap()
-	for i := uint64(0); i < 30; i++ {
-		bma.Add(i)
-		bmr.Add(i)
-	}
-	bmr.Optimize()
-	bmi := bma.Intersect(bmr)
-	err := bitmapsEqual(bmi, bma)
-	if err != nil {
-		t.Fatalf("expected intersection to equal array")
-	}
-	err = bitmapsEqual(bmi, bmr)
-	if err != nil {
-		t.Fatalf("expected intersection to equal run")
-	}
-}
-*/
 func TestShiftArray(t *testing.T) {
 	tests := []struct {
 		array []uint16
@@ -3883,6 +4079,47 @@ func BenchmarkUnionInPlaceRegression(b *testing.B) {
 	})
 }
 
+func BenchmarkUnionRunRunInPlace(bm *testing.B) {
+	bm.Skip("Skipping long running BenchmarkUnionRunRunInPlace")
+
+	runs := []struct {
+		name string
+		fn   func() []interval16
+	}{
+		{"FirstBitSet", runFirstBitSet},
+		{"LastBitSet", runLastBitSet},
+		{"FirstBitUnset", runFirstBitUnset},
+		{"LastBitUnset", runLastBitUnset},
+		{"InnerBitsSet", runInnerBitsSet},
+		{"OuterBitsSet", runOuterBitsSet},
+		{"OddBitsSet", runOddBitsSet},
+		{"EvenBitsSet", runEvenBitsSet},
+	}
+
+	for _, ar := range runs {
+		for _, br := range runs {
+			bm.Run("RunToBitmapRun-"+ar.name+"_"+br.name, func(bm *testing.B) {
+				for i := 0; i < bm.N; i++ {
+					arun := doContainer(containerRun, ar.fn())
+					brun := doContainer(containerRun, br.fn())
+
+					abmp := arun.runToBitmap()
+					unionBitmapRunInPlace(abmp, brun)
+				}
+			})
+
+			bm.Run("RunRun-"+ar.name+"_"+br.name, func(bm *testing.B) {
+				for i := 0; i < bm.N; i++ {
+					arun := doContainer(containerRun, ar.fn())
+					brun := doContainer(containerRun, br.fn())
+
+					unionRunRunInPlace(arun, brun)
+				}
+			})
+		}
+	}
+}
+
 func TestBitmapAny(t *testing.T) {
 	bm := NewBTreeBitmap()
 	if bm.Any() {