Here is a Python script that gives an upper bound of ~10^51.431 positions (legal and illegal). Moreover, the script doesn't check for legality of positions or if a position is reachable, so the proportion of legal positions should be a few orders of magnitude lower. A bootstrap approach can be used to approximate this proportion.
Note the following variables are updated:
squaresAvailable = Number of squares available to use
whiteProm = Maximum number in promotions available on the position for white pieces
blackProm = Maximum number in promotions available on the position for black pieces
import math
sum = 0
tosum = 43
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)
whiteProm = 8 - whitepawn
blackProm = 8
for blackpawn in range(9):
bp = math.comb(48 - whitepawn, blackpawn)
whiteProm = 8 - whitepawn
blackProm = 8 - blackpawn
squaresAvailable = 62 - whitepawn - blackpawn
for whitelsbishop in range(2 + max(whiteProm, 0)):
squaresAvailable = 62 - whitepawn - blackpawn - whitelsbishop
wbls = math.comb(math.ceil(squaresAvailable/2), whitelsbishop)
whiteProm = 8 - whitepawn - max(0, whitelsbishop - 1)
blackProm = 8 - blackpawn
for blacklsbishop in range(2 + max(blackProm, 0)):
squaresAvailable = 62 - whitepawn - blackpawn - whitelsbishop - blacklsbishop
bbls = math.comb(math.ceil(squaresAvailable/2), blackbishop)
whiteProm = 8 - whitepawn - max(0, whitelsbishop - 1)
blackProm = 8 - blackpawn - max(0, blacklsbishop - 1)
for whitedsbishop in range(2 + max(whiteProm, 0)):
squaresAvailable = 62 - whitepawn - blackpawn - whitelsbishop - blacklsbishop - whitedsbishop
wbds = math.comb(math.ceil(squaresAvailable/2), whitedsbishop)
whiteProm = 8 - whitepawn - max(0, whitelsbishop - 1) - max(0, whitedsbishop - 1)
blackProm = 8 - blackpawn - max(0, blacklsbishop - 1)
for blackdsbishop in range(2 + max(whiteProm, 0)):
squaresAvailable = 62 - whitepawn - blackpawn - whitelsbishop - blacklsbishop - whitedsbishop - blackdsbishop
bbds = math.comb(math.ceil(squaresAvailable/2), whitedsbishop)
whiteProm = 8 - whitepawn - max(0, whitelsbishop - 1) - max(0, whitedsbishop - 1)
blackProm = 8 - blackpawn - max(0, blacklsbishop - 1) - max(0, blackdsbishop - 1)
for whitenight in range(3 + max(whiteProm, 0)):
squaresAvailable = 62 - whitepawn - blackpawn - whitelsbishop - blacklsbishop - whitedsbishop - blackdsbishop - whitenight
wn = math.comb(squaresAvailable, whitenight)
whiteProm = 8 - whitepawn - max(0, whitelsbishop - 1) - max(0, whitedsbishop - 1) - max(0, whitenight - 2)
blackProm = 8 - blackpawn - max(0, blacklsbishop - 1) - max(0, blackdsbishop - 1)
for blacknight in range(3 + max(blackProm, 0)):
squaresAvailable = 62 - whitepawn - blackpawn - whitelsbishop - blacklsbishop - whitedsbishop - blackdsbishop - whitenight - blacknight
bn = math.comb(squaresAvailable, blacknight)
whiteProm = 8 - whitepawn - max(0, whitelsbishop - 1) - max(0, whitedsbishop - 1) - max(0, whitenight - 2)
blackProm = 8 - blackpawn - max(0, blacklsbishop - 1) - max(0, blackdsbishop - 1) - max(0, blacknight - 2)
for whiterook in range(3 + max(whiteProm, 0)):
squaresAvailable = 62 - whitepawn - blackpawn - whitelsbishop - blacklsbishop - whitedsbishop - blackdsbishop - whitenight - blacknight - whiterook
wr = math.comb(squaresAvailable, whiterook)
whiteProm = 8 - whitepawn - max(0, whitelsbishop - 1) - max(0, whitedsbishop - 1) - max(0, whitenight - 2) - max(0, whiterook - 2)
blackProm = 8 - blackpawn - max(0, blacklsbishop - 1) - max(0, blackdsbishop - 1) - max(0, blacknight - 2)
for blackrook in range(3 + max(blackProm, 0)):
squaresAvailable = 62 - whitepawn - blackpawn - whitelsbishop - blacklsbishop - whitedsbishop - blackdsbishop - whitenight - blacknight - whiterook - blackrook
br = math.comb(squaresAvailable, blackrook)
whiteProm = 8 - whitepawn - max(0, whitelsbishop - 1) - max(0, whitedsbishop - 1) - max(0, whitenight - 2) - max(0, whiterook - 2)
blackProm = 8 - blackpawn - max(0, blacklsbishop - 1) - max(0, blackdsbishop - 1) - max(0, blacknight - 2) - max(0, blackrook - 2)
for whitequeen in range(2 + max(whiteProm, 0)):
squaresAvailable = 62 - whitepawn - blackpawn - whitelsbishop - blacklsbishop - whitedsbishop - blackdsbishop - whitenight - blacknight - whiterook - blackrook - whitequeen
wq = math.comb(squaresAvailable, whitequeen)
whiteProm = 8 - whitepawn - max(0, whitelsbishop - 1) - max(0, whitedsbishop - 1) - max(0, whitenight - 2) - max(0, whiterook - 2) - max(0, whitequeen -1)
blackProm = 8 - blackpawn - max(0, blacklsbishop - 1) - max(0, blackdsbishop - 1) - max(0, blacknight - 2) - max(0, blackrook - 2)
for blackqueen in range(2 + max(blackProm, 0)):
squaresAvailable = 62 - whitepawn - blackpawn - whitelsbishop - blacklsbishop - whitedsbishop - blackdsbishop - whitenight - blacknight - whiterook - blackrook - whitequeen - blackqueen
bq = math.comb(squaresAvailable, blackqueen)
blackProm = 8 - blackpawn - max(0, blacklsbishop - 1) - max(0, blackdsbishop - 1) - max(0, blacknight - 2) - max(0, blackrook - 2) - max(0, blackqueen -1 )
sum = sum + wp * bp * wbls * bbls * bbds * wbds * 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: 51.43108306627925