Tuesday, March 11, 2014

GNU MP arbitrary precision library.

Tonight I am learning the arbitrary precision math library from gnu called gmp.   It is tough to work only with functions instead of just using math expressions.   It takes multiple lines of obscure code to do what you can do with a single C expression.  Instead of saying:
bb = floor (n *sqrt(.5) +1);
You have to say:
  mpf_t n;  mpf_init2(n, 256);  mpf_set_ui(n, 10);  mpf_pow_ui(n,n,12);
  mpf_t bb;  mpf_init2(bb, 256);  mpf_set_ui(bb, 5);  mpf_div_ui(bb, bb, 10);
  mpf_t srpf;  mpf_init2(srpf, 256);  mpf_sqrt(srpf, bb);

  mpf_mul(bb, srpf, n);  mpf_add_ui(bb,bb, 1);  mpf_floor(bb, bb);

Why deal with this extra complexity?  Because it is amazingly easy to run into problems of precision with just the 64 bits that the standard data structures give you.  Simple things like working with numbers in the 10^16 range, or 100!, or the 5000th Fibonacci number and you get the wrong values coming out of your functions.   The 5000th Fibonacci number runs you out of bits in the exponent of long double.

I am playing around with the problems on the Project Euler site and am running into a lot of situations where I cannot solve the problems with C because the problem sets have large enough numbers that they are even overflowing 64 bit integers and floating point numbers are fuzzy enough that two things that should match don't because the lowest precision bits don't match, when they should.

No comments:

Post a Comment