Make very small, (or large), exponential calculations

2
votes

Exponential limit of most 32 bit machines is 

exp( +/-700 )

But I would like to do an exponential calculation 

res = exp( x ) / exp( d )

when x or d are bigger than 700 I use the fact that 

exp( x + y ) = exp( x ) . exp( y )

So my calculation would be something along the line of

res = (exp( x - z ).exp(z)) / (exp( d - z ).exp(z))

or 

res = exp( x - z ) / exp( d - z )

where (x-z) < 700

But this approach is flawed in some cases, for example where x = 6000 and d = 10000
If we use z=5300 then

res = exp( 6000 - 5300 ) / exp( 10000 - 5300 )

res = exp( 700 ) / exp( 47000 )

But exp( 47000 ) = 0 on a 32 bit machine.

If I replace z = 9300 then I get the opposite effect.

res = exp( -3300 ) / exp( 700 )

So how could I solve the above equations, (that should return a 32bit valid number I think), given the limitations of the computer?

Edit The reason for doing this is I am using the formula

P( a ) = P(y1) * P(y2) * P(y3) ... P(yN)

In order to prevent overflow I then do

a = log( P(y1) ) + log( P(y2) ) +  log (P(y3)) ... log( P(yN) )

b = log( P(z1) ) + log( P(z2) ) +  log (P(z3)) ... log( P(zN) )

...

z = log( P(zz1) ) + log( P(zz2) ) +  log (P(zz3)) ... log( P(zzN) )

to get a total I do 

total = a + b ... z

and to calculate the percentage I do 

(exp(a) / exp( total ) ) * 100

but it is possible that "a" and/or "total" are greater than 700

I guess the question could be how could I calculate the percentage without using the exponential

3
c++math32-bitexponential
you do realize that exp(x-y) = exp(x)/exp(y) which would mean that exp(6000)/exp(10000) = exp(-4000), which is past the limitations – R Nar
I cannot really follow you math. exp(x) / exp(d) with x=6000 and d=10000 is exp(6000-10000)=exp(-4000) and this wont produce overflow – 463035818_is_not_a_number
Search internet for "Big Number Library C++". – Thomas Matthews
What is the context of this computation? Usually, when extremely large numbers appear in the numerator and denominator of a fraction, but the fraction itself is readily representable, it is time to re-order the math. If you think the math cannot be re-arranged, you could add an additional scale factor to each of your operands, creating a custom "uber-float" format with two exponents. – njuffa
I added some more information to my question to hopefully give more details. – Simon Goodman

3 Answers

4
votes

It doesn't matter that the answer should be a 32 bit number if some of the intermediate steps in the calculations aren't.

For math that goes outside the bounds of an int or long type, you probably need to start using something like GMP.

https://gmplib.org/

3
votes

I assume that you want to compute this:

p = exp(a) / exp(b)

And since a^b/a^c == a^(b-c) this reduces to

p = exp(a - b)

which can be easily computed if that difference is below that critical exponent.

If it isn't, then your result cannot be represented by primitive datatypes like double (because it's either extremely large or extremely small), you then need some kind of arbitrary precision numbers, probably provided by some library.

But if you only need to print the result, or store it somehow, then you can easily compute even extremely large numbers:

For that, you change to base 10 (for displaying), compute the equivalent exponent therefore (tExponent = log10(eExponent)), and get that value into the allowed range between std::numeric_limits::max_exponent10 and std::numeric_limits::min_exponent10, saving the difference as scaling factor.

For now, I just have a quick and dirty live example, showing

exp(90000) / exp(100)  =  1.18556 scaled by 10^39043

(Check at wolfram alpha)

Note: When I wrote this, it was pretty late in the evening. I'm leaving this here for an "alternative" approach.

Now, generally, there's

a^b = [a^(b/c)]^c

And since

(a/b)^c = (a^c)/(b^c)

holds, too, I guess the easiest approach here is to just divide both exponents as long as one of them is above your critical value, then do the exponentiation, divide the results, and finally use the divisor of the former exponents as exponent for the quotient:

double large_exp_quot(
    double eNum,
    double eDenom,
    unsigned int const critical = 200) {
  if (abs(eNum - eDenom) > critical) {
    throw out_of_range{"That won't work, resulting exponent is too large"};
  }
  unsigned int eDivisor = 1;
  while (abs(eNum) > critical or abs(eDenom) > critical) {
    eNum /= 2;
    eDenom /= 2;
    eDivisor *= 2;
  }
  return pow(exp(eNum) / exp(eDenom), eDivisor);
}

But this will only work, if the result of your computation can actually be stored using the C++ primitive datatypes, in this case double. The example you gave ... with exponents 6000 and 10000 ... is obviously not representable with a double (it's e^(-4000) and thus incredibly small)

2
votes

Numerically unstable computation: exp(a) / exp(b)
Equivalent stable computation: exp(a - b)

Numerically unstable computation: Πi=1..n pi
Equivalent stable computation: exp(Σi=1..n log(pi))