This approach:
train it as a classification problem which output confidence, and hope that the confidence will reflect the actual probability. Sometimes the network would output the correct probability and fail the test, but on average it would tend to the correct probability.
will work with some limitations. If you use a classifier with softmax activation and multiclass log loss
$$\mathcal{L}(\mathbf{\hat{y}},\mathbf{y}) = -\mathbf{y} \cdot \text{log}(\mathbf{\hat{y}})$$
where $\mathbf{\hat{y}}$ is the network output as a vector, and $\mathbf{y}$ is the actual output from an individual sample. Your input should be the settings of the die.
Optimising this loss will converge on approximated probabilities for each discrete output. You can demonstrate this with some simple examples - for instance if you train a network with a single input - a one-hot-encoded die type from the classic D&D dice sets, plus deliberately chosen examples of different results in the right frequencies, you will end up with a classifier that predicts roughly $p=0.25$ for results of 1,2,3,4 for a d4 and $p=0.125$ for results of 1,2,3,4,5,6,7,8 for a d8
So it works mathematically. Whether it works for your situation depends on details. You need enough data samples to cover both the distribution of results under each setting, and any complexities of how that distribution varies with the settings. In the limit of wanting very accurate predictions of probability within a complex space you will need a huge number of samples. You should be able to find a compromise between accuracy and generalisation by trying different levels of regularisation - this will be necessary as over-fitting to input/output sample pairs as seen is going to be a serious problem for a neural network trained on this data.
One thing you can do to help a classifier learn probabilities is always take some number of samples with the same settings - e.g. 10 or 100 or 1000 each with the same settings - that should guarantee that the network cannot simply converge to predict high $p$ values for single outputs as seen, as it will have counter-examples to work with.
You mention that you have 40 dimensions of settings. Whether this is an issue will depend on how the probability distribution varies based on those settings. However, at minimum you should be thinking in terms of millions of samples for training, or possibly a fast on-demand generator that can generate 1000s of new samples per second to train with.
You can test accuracy by building histograms using some fixed (and as yet unseen) setting and comparing to NN predictions of probabilities of that setting. Even getting accurate test results is likely to require 1000s of samples.
However it may require a lot of training and data.
If you cannot obtain a very large training set here, then a purely statistical "black box" approach is probably not feasible, regardless of whether you use neural networks, or more raw analysis. What neural networks add is smooth interpolation between different settings values, as a form of approximation. This seems desirable for your problem, as you will never fully explore 40 dimensions of 100 values in the lifetime of the universe - but you need some confidence that minor changes in settings equate to minor changes in probability distributions in most cases.
It's OK to have one or two major shifts across the input space, but if the distribution depends in some cryptographic primitive or similar complex high frequency (over space) and high amplitude on the input variables, there is no way to obtain approximations using statistics.
The alternative to statistical approaches is to find some way to break open this black box through analysis. No AI system can do that in general at the moment, so you would rely on human ingenuity.