By "bm rng" I mean the Method used to calculate the random number. It's called "Box–Muller transform", I just forgot the name and shortened it.
I called stddev that is because it is the same calculation. It is halve the rangesize and thus the distance between mean and min/max of the range.
The way the final stat value is calculated is weird or I might be misunderstanding something. Both stddev and mean are divided by 1000 and thus normalized to the max range (0=0 ... 1000=1). rir_bm() calls random_bm() with stddev and mean and returns a factor which is then multiplied by the rangesize and adds that to the start value. The problem is the multiplication with the rangesize.
If you look at the distribution of the values random_bm() returns, you get a normal distribution. The thing is that all these values are normalized to the max range of 1000. If you simply multiply them by 1000 you get a normal distribution of stat values centered correctly around mean, which I thought was what we are trying to calculate here.
The important part is that the the values random_bm() returns are absolute values(with a root of 0), but rangesize is relative. It's a bit difficult for me to explain, but basically this rangesize can be anywhere on the range. Similar to your example, if you take values on the higher end 800 to 900 for example. The mean is 850, but if you calculate it with the code you get num=0(mean of the box-muller transform)*stddev+mean=0.85; diff=rangesize*num=100*0.85=85; value=start+diff=800+85=885. And well 850=/=885.
And since both rir_bm and random_bm are basically trying to calculate the same thing in two different ways, I thought just make it simple and use the method I described.