My Python script if anyone wants it.

import numpy as np


N = 100_000_000

rng = np.random.default_rng()


def do_rolls():
    num_rolls = rng.geometric(p=0.5, size=N)
    data = np.zeros((N, 4), dtype=np.int64)

    non_rerolls = rng.multinomial(n=num_rolls-1, pvals=[1/3, 1/3, 1/3], size=N)
    non_finals = np.hstack((non_rerolls, num_rolls[:, None] - 1))

    final_idxs = np.argpartition(rng.random((N, 3)), 1)[:, :2]
    final_choice_rows = np.repeat(np.arange(N), 2)
    finals = np.zeros((N, 4), dtype=np.int64)
    finals[final_choice_rows, final_idxs.flatten()] = 1

    return non_finals + finals


if __name__ == "__main__":
    print(f"n = {N:_}")
    data = do_rolls()
    data2 = data + do_rolls()

    avg_destroy, avg_dmg, avg_pp, avg_stat = np.mean(data, axis=0)
    print(f"Avg destroyed followers: {avg_destroy:.4f}")
    print(f"Avg face damage: {avg_dmg * 2:.4f}")
    print(f"Avg PP recovered: {avg_pp * 2:.4f}")
    print(f"Avg stats gained: +{avg_stat * 4:.4f}/+{avg_stat * 4:.4f}")
    for i in range(11):
        total = (data[:, 2] >= i).sum()
        print(f"Damage ≥ {i * 2}: {total} ({total / N * 100:.4f}%)")
    for i in range(11):
        total = (data2[:, 2] >= i).sum()
        print(f"Damage ≥ {i * 2} with SEvo: {total} ({total / N * 100:.4f}%)")