开发者

c++ log function is using floating point precision

I am having an interesting seg fault in the following function when I give it a number very close to 1.0. Specifically when the number would be rounded to 1.0 at FLOATING POINT precision.

double get_random_element(double random_number)
{
    if (random_number <= 0.0 || random_number >= 1.0)
        throw std::runtime_error("Can't have a random number not on the range (0.0, 1.0)");
    return -log(-log(random_number));
}

If random_number is 1.0 then log(1.0) = 0.0 and the log of zero is an undefined calculation leading to a seg fault. However I would have thought that the error checking on the first line would have prevented this from ever happening. Ddebugging shows that a number very close to 1 will pass through the error checking but return 0 from the log function anyway leading me to believe that the log function is using only single floating point precision.

my includes are as follows so i can only assume I'm using the log from math.h

#include <string>
#include <math.h>
#include <sstream>
#include <map>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/uniform_int.hpp>
#include <boost/random/uniform_real.hpp>
#include <boost/random/variate_generator.hpp>
#include <utility>

UPDATE: As pointed out an easy solution is to just use a floating point number as the argument and if a number equal to 1.0f is passed in to just remove std::numeric_limits::epsilon() to give a number which can be safely passed into the double log.

But the question I'd like answered is why does calling double log of a number near but not equal to 1 fail.

UPDATE 2: After recreating this problem in a test project I think the problem is actually in the inputs. If I pass in

double val = 1.0 - s开发者_开发问答td::numerical_limits<double>::epsilon();

I have no problems with the function. However what I actually pass in is

boost::mt19937 random_generator;
double val = (random_generator()+1)/4294967297.0;

where random_generator is designed to return a number on the range [0, 2^32 - 1] == [0,4294967295]. So I decided to punch in the largest possible return value

double val = (4294967295+1)/4294967297.0;

which quickly gave me a warning about unsigned int overflow and sure enough generated a zero. I am recompiling with the following:

get_random_element((random_generator()+1.0)/4294967297.0);

and hopefully this strange behaviour will be resolved.

UPDATE 3: I have finally found what is going on here... and as usual it comes down to user error (myself being the error). There was a second control path leading to this method which temporarily stored the double value as a float and then converted it back to double leading to 0.999999999 being rounded to 1.0 and then passed into the -log(-log(x)) function and causing it to fall over. What I still don't understand is why my checking

 if (random_number <= 0.0 || random_number >= 1.0) throw runtime_error(blah)

didn't catch the erroneous input before it was passed into the log functions?


I think quamrana has a good point (it immediately drew my attention too). However, I've been able to run this snippet for considerable length:

#include <math.h>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/uniform_real.hpp>

double get_random_element(double random_number)
{
    if (random_number <= 0 || random_number >= 1.0f)
        throw std::runtime_error("Can't have a random number not on the range (0.0, 1.0)");
    return -::log(-::log(random_number));
}

int main()
{
    boost::mt19937 rng; 
    boost::uniform_real<> random(std::numeric_limits<double>::epsilon(),1);
    for (;;)
    {
        double r = random(rng);
        double gre = get_random_element(r);
        std::cout << "r = " << r << ", gre = " << gre << std::endl;
    }
    return 0; // not reached
}

E.g.:

sehe@meerkat:/tmp$ ./t | grep '^r = 0.999999' 
r = 0.999999, gre = 14.4777
r = 0.999999, gre = 13.7012
r = 0.999999, gre = 14.0492
r = 0.999999, gre = 14.1161
[.... many many lines snipped ....]
r = 0.999999, gre = 14.3691
r = 0.999999, gre = 13.424
r = 0.999999, gre = 14.4822
r = 0.999999, gre = 14.286
r = 0.999999, gre = 14.4344
r = 0.999999, gre = 14.0572
r = 0.999999, gre = 14.0607
r = 0.999999, gre = 14.1126
r = 0.999999, gre = 13.575
r = 0.999999, gre = 13.4754
r = 0.999999, gre = 13.5486
r = 0.999999, gre = 14.1983
^C

real    18m14.005s
user    20m5.667s
sys 12m19.302s

Perhaps you could use something similar in vein?

0

上一篇:

下一篇:

精彩评论

暂无评论...
验证码 换一张
取 消

最新问答

问答排行榜