calculating exp(-300)

[Nathan Moore]

Hello list,

I'm working on a monte-carlo simulation which requires the calculation 
of probability for rather high energy levels, Boltzmann probability 
defined, exp(-u/T).  The quantity u/T is occasionally very large, and 
in some cases always large.  I'd like to use GMP rational numbers 
(mpq_t) to represent the probability and partition function of the 
system.

While the standard function, exp(), doesn't seem to exist for rational 
numbers, I see that an alternative form, mpq_div_2exp() does exist.  I 
feel like there should be some way to represent exp(-a) as y / 2^x, but 
with the constraint that x be an integer the best method of 
representation isn't immediately obvious to me.

At present I can write down that,
	-a = -x ln(2) + ln(y).
and then estimate x by ignoring y,
	x = a / ln(2)
and then round down x, and find the correction y supplies,
	x = floor(x);
	y = exp(-a + x*ln(2)).
Then with base 2 representation I can specify the weight (y converted 
to mpq_t),
	mpq_div_2exp(weight, y, x)

I think this will work, but I worry it is too primitive.  Has anyone on 
the list addressed this problem differently?

regards,

Nathan Moore
University of Minnesota, Physics