Dice rolls stats

def v1():
    "Your original version"

    i = j = k = l = m = n = 1
    total = maximumRolls = 0
    while i < 7:
        while j < 7:
            while k < 7:
                while l < 7:
                    values = [i, j, k, l]
                    values.sort()
                    subtotal = values[1] + values[2] + values[3]
                    total += subtotal
                    if subtotal == 18:
                        maximumRolls += 1
                    l += 1
                k += 1
                l = 1
            j += 1
            k = 1
        i += 1
        j = 1

    total = total / 1296

    print("Average of 4d6kh3 = " + str(total)[:5] + " Max Rolls: " + str(maximumRolls) + "/1296 (" + str(maximumRolls * 100 / 1296)[:4] + "%)")



from itertools import product
def v2(num_dice, faces_per_die, num_highest_to_keep):
    "itertools.product to generate all the possible rolls"

    max_roll_count = total_sum = 0
    total_combinations = faces_per_die ** num_dice

    for roll in product(range(1, faces_per_die + 1), repeat=num_dice):
        sorted_roll = sorted(roll, reverse=True)
        kept_sum = sum(sorted_roll[:num_highest_to_keep])
        total_sum += kept_sum

        if kept_sum == num_highest_to_keep * faces_per_die:
            max_roll_count += 1

    average = total_sum / total_combinations
    percentage = (max_roll_count / total_combinations) * 100

    print(f"Average of {num_dice}d{faces_per_die}kh{num_highest_to_keep} = {average:.2f} "
          f"Max Rolls: {max_roll_count}/{total_combinations} ({percentage:.2f}%)")

import numpy as np

def v3(num_dice, faces_per_die, num_highest_to_keep):
    """Recursively generating the the rolls, we start at the top with no roll
    and we go further adding one roll value in the roll at a time,
    when n==0 we are at the bottom layer of the recursion and we have our complete roll with n num_dice values
    and we propagate the 2 values we care about to the top : total_sum and max_roll_count
    by doing the sum of all the n-1 layers values until we are back to the top
    """
    def v3_rec(n, roll):
        if n == 0:
            sorted_roll = sorted(roll, reverse=True)
            kept_sum = sum(sorted_roll[:num_highest_to_keep])
            max_roll = 1 if kept_sum == num_highest_to_keep * faces_per_die else 0
            return [kept_sum, max_roll]

        return np.array([v3_rec(n-1, roll+[i]) for i in range(1, faces_per_die+1)]).sum(axis=0)

    total_sum, max_roll_count = v3_rec(num_dice, [])
    total_combinations = faces_per_die ** num_dice

    average = total_sum / total_combinations
    percentage = (max_roll_count / total_combinations) * 100

    print(f"Average of {num_dice}d{faces_per_die}kh{num_highest_to_keep} = {average:.2f} "
          f"Max Rolls: {max_roll_count}/{total_combinations} ({percentage:.2f}%)")


def v4(num_dice, faces_per_die, num_highest_to_keep):
    """From a roll we generate a new roll, for example with 4d6, if we are at [1,1,1,1] the
    next roll is [1,1,1,2], and if we are at [1,1,1,6] the next one is [1,1,2,1]
    and once we are at [6,6,6,6] then we know its finished
    """
    max_roll_count = total_sum = 0
    total_combinations = faces_per_die ** num_dice

    roll = [1 for _ in range(num_dice)]
    returned_to_start = False
    while not returned_to_start:
        sorted_roll = sorted(roll, reverse=True)
        kept_sum = sum(sorted_roll[:num_highest_to_keep])
        total_sum += kept_sum
        if kept_sum == num_highest_to_keep * faces_per_die:
            max_roll_count += 1

        carry_the_addition = True
        idx = num_dice-1

        while carry_the_addition:
            if roll[idx] + 1 > faces_per_die:
                if idx == 0:
                    returned_to_start = True
                    break

                roll[idx] = 1
                idx -= 1
            else:
                roll[idx] += 1
                carry_the_addition = False

    average = total_sum / total_combinations
    percentage = (max_roll_count / total_combinations) * 100

    print(f"Average of {num_dice}d{faces_per_die}kh{num_highest_to_keep} = {average:.2f} "
          f"Max Rolls: {max_roll_count}/{total_combinations} ({percentage:.2f}%)")

def v5(num_dice, faces_per_die, num_highest_to_keep):
    """Same as v4 but with for example 4d6, the rolls are from 0 to 5 when doing
    next roll calculations so we can use % and // operators more easily,
    to get the true roll we just add 1 to every value of the roll
    """

    max_roll_count = total_sum = 0
    total_combinations = faces_per_die ** num_dice
    for i in range(total_combinations):
        roll_minus_one = []
        temp = i
        for _ in range(num_dice):
            roll_minus_one.append(temp % faces_per_die)
            temp //= faces_per_die

        roll = [v+1 for v in roll_minus_one]
        sorted_roll = sorted(roll, reverse=True)

        kept_sum = sum(sorted_roll[:num_highest_to_keep])
        total_sum += kept_sum
        if kept_sum == num_highest_to_keep * faces_per_die:
            max_roll_count += 1

    average = total_sum / total_combinations
    percentage = (max_roll_count / total_combinations) * 100
    print(f"Average of {num_dice}d{faces_per_die}kh{num_highest_to_keep} = {average:.2f} "
          f"Max Rolls: {max_roll_count}/{total_combinations} ({percentage:.2f}%)")


def v6(num_dice, faces_per_die):
    """I asked the LLM o4-mini to come up with a mathematical formula to compute things
    i asked the formula for the special case when we take the highest n-1 die because its easier to compute
    Because in this case the expectancy of the n-1 highest dice is the expetancy of the n dice
    minus the expectancy of the min value of the n dice.
    Function generated by o4-mini
    """
    num_highest_to_keep = num_dice - 1
    n, f = num_dice, faces_per_die

    # 1) E[min of n dice] = sum_{k=1}^f P(min >= k)
    #                  = sum_{k=1}^f ((f+1-k)/f)^n
    E_min = sum(((f + 1 - k) / f) ** n for k in range(1, f + 1))

    # 2) E[sum of all n dice] = n * (f+1)/2
    #    so E[sum of top n-1] = n*(f+1)/2 - E_min
    average = n * (f + 1) / 2 - E_min

    # 3) P(max roll) occurs if at least n-1 dice = f
    #    = C(n,n-1)*(1/f)^(n-1)*((f-1)/f) + (1/f)^n
    #    = (n*(f-1) + 1) / f^n
    p_max = (n * (f - 1) + 1) / (f**n)

    print(f"Expectancy of {num_dice}d{faces_per_die}kh{num_highest_to_keep} = {average:.2f} "
          f"Max Rolls expectancy : {p_max*100:.2f}%")




nb_dices = 4
nb_dice_faces = 6
top_n = 3

v1()
v2(nb_dices, nb_dice_faces, top_n)
v3(nb_dices, nb_dice_faces, top_n)
v4(nb_dices, nb_dice_faces, top_n)
v5(nb_dices, nb_dice_faces, top_n)
v6(nb_dices, nb_dice_faces)