Incorrect upper bound on number of possible chess positions (10^41)?

Question

I've written a script implementing a basic idea to count the number of chess positions. I've not even accounted for legality of the positions or difference in light and dark square bishops. Despite this, my script only counted ~10^40.7 possible chess positions. Accounting for the previous discrepancies, this program can easily reduce the total positions to 10^38, which is an underestimation compared to the upper bound of 8.7 * 10^45.

I keep track of three general variables:

squaresAvailable = Number of squares available to use
whiteProm = Maximum number of promotions available on the position for white pieces
blackProm = Maximum number of promotions available on the position for black pieces

I iterate over each possibility of number of white and black pawns, bishops, rooks, and queens. For example, if there three white queens on the board, the number of possibilities that those queens can specifically be arranged is wq = (squaresAvailable choose 3).

Here is the script (Python):

import math
sum = 0

for whitepawn in range(9): #whitepawn is the number of white pawns on the board (0 to 8)
    #wp represents the number of ways to arrange a given number of white pawns
    wp = math.comb(48, whitepawn)
    
    for blackpawn in range(9):
        bp = math.comb(48 - whitepawn, blackpawn)
        whiteProm = 8 - whitepawn
        blackProm = 8 - blackpawn
        squaresAvailable = 62 - whitepawn - blackpawn
        
        for whitebishop in range(3 + max(whiteProm, 0)):
            wb = math.comb(squaresUsed, whitebishop)
            squaresAvailable =  62 - whitepawn - blackpawn - whitebishop
            
            if whitebishop > 2:
                #if there are more than 2 white bishops, some of them must've come from 
                #promotion, reducing the rest of the available white promotions
                
                whiteProm = whiteProm - (whitebishop - 2)
            
            for blackbishop in range(3 + max(blackProm, 0)):
                bb = math.comb(squaresUsed, blackbishop)
                squaresAvailable = 62 - whitepawn - blackpawn - whitebishop - blackbishop
                if blackbishop > 2:
                    blackProm = blackProm - (blackbishop - 2)
                
                for whitenight in range(3 + max(whiteProm, 0)):
                    wn = math.comb(squaresUsed, whitenight)
                    squaresAvailable = 62 - whitepawn - blackpawn - whitebishop - blackbishop - whitenight
                    if whitenight > 2:
                        whiteProm = whiteProm - (whitenight - 2)

                    for blacknight in range(3 + max(blackProm, 0)):
                        bn = math.comb(squaresUsed, blacknight)
                        squaresAvailable = 62 - whitepawn - blackpawn - whitebishop - blackbishop - whitenight - blacknight
                        if blacknight > 2:
                            blackProm = blackProm - (blacknight - 2)

                        for whiterook in range(3 + max(whiteProm, 0)):
                            wr = math.comb(squaresUsed, whiterook)
                            squaresAvailable = 62 - whitepawn - blackpawn - whitebishop - blackbishop - whitenight - blacknight - whiterook
                            if whiterook > 2:
                                whiteProm = whiteProm - (whiterook - 2)
                            
                            for blackrook in range(3 + max(blackProm, 0)):
                                br = math.comb(squaresUsed, blackrook)
                                squaresAvailable = 62 - whitepawn - blackpawn - whitebishop - blackbishop - whitenight - blacknight - whiterook - blackrook
                                if blackrook > 2:
                                    blackProm = blackProm - (blackrook - 2)

                                for whitequeen in range(2 + max(whiteProm, 0)):
                                    wq = math.comb(squaresUsed, whitequeen)
                                    squaresAvailable = 62 - whitepawn - blackpawn - whitebishop - blackbishop - whitenight - blacknight - whiterook - blackrook - whitequeen
                                    if whitequeen > 1:
                                        whiteProm = whiteProm - (whitequeen - 1)

                                    for blackqueen in range(2 + max(blackProm, 0)):
                                        bq = math.comb(squaresUsed, blackqueen)
                                        squaresAvailable = 62 - whitepawn - blackpawn - whitebishop - blackbishop - whitenight - blacknight - whiterook - blackrook - whitequeen - blackqueen
                                        if blackqueen > 1:
                                            blackProm = blackProm - (blackqueen - 1)
                                    
                                        sum = sum + wp * bp * wb * bb * wn * bn * wr * br * wq * bq
                                    
#64 squares available for white king and 63 for black king as an upper bound
print(math.log10(sum * 64 * 63))

OUTPUT: 40.73915095995693

As previously mentioned, I think I can lower this number by a significant amount considering for legality, dark square and light square bishops, so there must've been an error in my code somewhere. Where did my code go wrong?

EDIT: I found the error, I wasn't updating whiteProm and blackProm properly after each loop. After fixing this issue, the number of chess positions returned is only 10^45.5, which still seems very low.

EDIT2: I've worked out the other errors, incorporated the difference between light and dark square bishops and my upper estimate on the number of chess positions is 10^51, so the program output is correct.