diff --git a/roaring/roaring.go b/roaring/roaring.go
index 24e94221b..0bd4da63a 100644
--- a/roaring/roaring.go
+++ b/roaring/roaring.go
@@ -3141,8 +3141,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)
 		}
 	}
 	panic(fmt.Errorf("invalid union op: unknown types %d/%d", c.typ(), other.typ()))
@@ -4712,6 +4711,198 @@ func unionBitmapBitmapInPlace(a, b *Container) *Container {
 	return a
 }
 
+// unions run b into run a, mutating a in place.
+func unionRunRunInPlace(a, b *Container) *Container {
+	statsHit("unionInPlace/RunRun")
+
+	a = a.Thaw()
+	runs, n := unionInterval16InPlace(a.runs(), b.runs())
+
+	a.setRuns(runs)
+	a.setN(n)
+	return a
+}
+
+// unionInterval16InPlace merges two slice of intervals 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
+
+		// subindex (ii) to state mapping
+		// .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 91168f656..085b6e1d1 100644
--- a/roaring/roaring_internal_test.go
+++ b/roaring/roaring_internal_test.go
@@ -721,6 +721,188 @@ func TestUnionMixed(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 TestIntersectMixed(t *testing.T) {
 	a := NewContainerRun([]interval16{{start: 5, last: 10}})
 	b := NewContainerArray([]uint16{1, 4, 5, 7, 10, 11, 12})
@@ -4149,3 +4331,77 @@ func TestDifferenceInPlace_N(t *testing.T) {
 		t.Error("expected difference of containers to have n=0")
 	}
 }
+
+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 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)
+				}
+			})
+		}
+	}
+}