Computing better average for the number of hits re: https://github.com/GameCrux/ManyWorlds and https://forum.maplelegends.com/index.php?threads/average-hits-calculator.38552

main.py

from constraint import Problem, MinSumConstraint
from functools import partial
from random import randint
from math import ceil


def solve_naive(lower, upper, hp):
    avg = (upper + lower) / 2
    return hp / avg


def _valid_min_hits(hp, *hits):
    """Each variable represents a hit on a monster, which can be within a range.
    We assert that the value is valid if we cannot drop one of the values and
    still maintain our original constraint."""
    acc = 0
    for i, hit in enumerate(hits):
        acc += hit
        if acc >= hp and i < len(hits) - 1:
            # the monster would have already been dead by this point
            return False
    return True


# return the number of solutions that are possible
def _solve_exact_k(lower, upper, hp, k):
    problem = Problem()
    vars = [f"hit_{i}" for i in range(k)]
    for var in vars:
        problem.addVariable(var, list(range(lower, upper + 1)))
    problem.addConstraint(MinSumConstraint(hp))
    problem.addConstraint(partial(_valid_min_hits, hp), vars)
    return problem.getSolutions()


def get_num_solutions(att_lower, att_upper, hp, bound_lower, bound_upper):
    # solution at index 0 is always 0
    res = [0] * bound_lower
    for i in range(bound_lower, bound_upper + 1):
        res_k = _solve_exact_k(att_lower, att_upper, hp, i)
        res.append(len(res_k))
    return res


def solve_exact(lower, upper, hp):
    min_hits = ceil(hp / upper)
    max_hits = ceil(hp / lower)
    # be aware, this will be very slow if the minimum number of hits to kill is large...
    problem = Problem()
    s = get_num_solutions(lower, upper, hp, min_hits, max_hits)
    print(s)
    # inclusive range
    r = upper - lower + 1

    # possible ways is based on the size of the range...
    return sum([i * s_i / (r ** i) for i, s_i in enumerate(s)])


def solve_approx(lower, upper, hp, samples=1000):
    min_hits = ceil(hp / upper)
    max_hits = ceil(hp / lower)

    # this may be slow if range of hits is wide
    weights = [0] * min_hits
    for k in range(min_hits, max_hits + 1):
        valid = 0
        for _ in range(samples):
            # sample uniformly from the range
            hits = [randint(lower, upper) for _ in range(k)]
            # check if this is valid
            valid += 1 if sum(hits) >= hp and _valid_min_hits(hp, *hits) else 0
        weights += [valid / samples]
    print(weights)
    return sum([i * w_i for i, w_i in enumerate(weights)]) / sum(weights)


def debug_exact(*args):
    print("exact", *args, "=>", solve_naive(*args), solve_exact(*args))


def debug_approx(*args):
    print("approx", *args, "=>", solve_naive(*args), solve_approx(*args))


debug_exact(40, 60, 150)
debug_approx(40, 60, 150)

debug_exact(40, 80, 150)
debug_approx(40, 80, 150)

# debug_exact(40, 120, 150)
debug_approx(1500, 3000, 15000)
debug_approx(1500, 5000, 15000)
debug_approx(1500, 8000, 15000)

requirements.txt

python-constraint==1.4.0

[0, 0, 0, 4796, 93765]
exact 40 60 150 => 3.0 3.48212935968038
[0, 0, 0, 0.517, 0.476]
approx 40 60 150 => 3.0 3.479355488418933
[0, 0, 66, 61255, 203360]
exact 40 80 150 => 2.5 3.03270411050333
[0, 0, 0.048, 0.889, 0.078]
approx 40 80 150 => 2.5 3.0295566502463047
[0, 0, 0, 0, 0, 0.0, 0.094, 0.671, 0.248, 0.008, 0.0]
approx 1500 3000 15000 => 6.666666666666667 7.166503428011754
[0, 0, 0, 0.0, 0.176, 0.533, 0.284, 0.039, 0.003, 0.0, 0.0]
approx 1500 5000 15000 => 4.615384615384615 5.188405797101449
[0, 0, 0.015, 0.404, 0.41, 0.124, 0.019, 0.0, 0.0, 0.0, 0.0]
approx 1500 8000 15000 => 3.1578947368421053 3.7201646090534983