Generate 128 bit Poisson distributed numbers in cpp

I try to measure probabilities in the order of 1e-25 or even lower. I run a cpp code which is basically an ODE solver numerically (with simple Euler steps). I have long long ints for storing the cell numbers (every step I roll a dice and I update the cell numbers using a 64 bit precision Mersenne Twister random num generator and I choose numbers from a poisson distribution with the mersenne engine). As long long int and 64 bit mersenne and poisson dist is only enough until approx. 1e19, I would need to have 128 bit precision. I am searching for a reliable method to generate 128 bit pseudo random numbers and also to store them like uint_128 and so, and generate especially 128bit random poisson distributed numbers.

I have searched in the topic but did not find anything helpful, maybe I could have puzzle the pieces together, but I could not...

Thank you in advance!