Are there any other (more efficient) methods (than the one I used) to solve my problem on Python?

So I was facing this math problem;

The essence of the problem goes as following: I am trying to calculate the amount of (whole) numbers in a given interval that can be expressed as the sum of multiples of given numbers.

In this case, running from 1 to 100, I want to find the numbers that can be written as 7x+13y, when x and y are whole numbers. (post scriptum: it has to do with donuts, hence the donuts mention in my code).

I am not an expert at Python, but I've managed to come up with something, by building a function that, for an interval of numbers [a,b] (in my case from 1 to 100) and specific numbers {n1,n2} (in my case 7 and 13), computes the numbers in the interval that can be written as n1 times j + n2 times k. But it takes 800 calculations to check for everything and I was wondering if there could be any other methods, more efficient ones, to solve for this kind of problem. Here is my code, annotated (I've put my final result at the end):

import numpy as n
import matplotlib.pyplot as plt
import os
from itertools import islice 
import math

#consider 13j+7k=i
#k represents the number of boxes of 7 donuts needed along j boxes of 13 donuts to obtain i donuts total
#i is between 1 and 100, j has to be < 8 since 8*13>100 but 7*13<100; so j is between 0 and 7
#we let i run between 1 and 100, j from 0 to 7, and for each combination we compute k
#we then order the list in sublists
#remove the negative values, keep whole ones
#make a list called results which compiles the number i of donuts for each solution found in k
#remove redundant solutions in a final_results list; print that list to display the solutions and check
#probabilities of fulfilling a random order of donuts chosen between 1 and 100 is the length of the final_results' list divided by a hundred (or b-a+1)

def function(a,b,n1,n2): #interval from a to b, with multiples n1 and n2
    k=[]
    results=[]
    j_range=math.floor((b-a+1)/n1)+1 #range for j
    
    for i in range((b-a)+2):
        if i>0:
            for j in range(j_range):
                k.append((i - n1*j)/n2) #here the k list is formed

    length_to_split=[j_range for i in range((b-a+1))]

    Input = iter(k) 
    Output = [list(islice(Input, elem)) for elem in length_to_split] #k list is ordered in sublists

    for i in range(len(Output)):
        for j in range(len(Output[i])):
            if float(Output[i][j]).is_integer() is True and Output[i][j]>=0:
                results.append(i+1) #results list is formed
                
    final_results = list(dict.fromkeys(results)) #remove possible redundant solutions
    print(final_results) #print the list of numbers that satisfy initial conditions

    probabilities = 100*(len(final_results)/(b-a+1))

    return print('probabilities of fulfilling the order is of', probabilities, '%')




function(1,100,13,7)

[7, 13, 14, 20, 21, 26, 27, 28, 33, 34, 35, 39, 40, 41, 42, 46, 47, 48, 49, 52, 53, 54, 55, 56, 59, 60, 61, 62, 63, 65, 66, 67, 68, 69, 70, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100]
probabilities of fulfilling the order is of 64.0 %

(I think my result is right, though I am not 100% sure it is correct).

So with the loops going from 1 to 100 and 0 to 7 the program has to do 800 calculations. I was wondering if there were any other more efficient methods to obtain the solution?

Thank you very much.

By the way my past question on this topic was closed despite the fact that I had edited the problems mentioned by the user who closed my question. That's not very urbane of you, kind sir.