Re: [boost] [review] Multiprecision review (June 8th - 17th, 2012)

Subject: Re: [boost] [review] Multiprecision review (June 8th - 17th, 2012)
From: Steven Watanabe (watanabesj_at_[hidden])
Date: 2012-06-25 00:33:08

AMDG

I really haven't had time to go through this
library as thoroughly as I would have liked,
but, since it's the last day of the review:
*I vote to accept Multiprecision into Boost*

Here are my comments on the code so far.
It's not very much, since I only got
through a few files.  I'll try to expand
on this as I have time, but no promises.

Revision: sandbox/big_number_at_79066

concepts/mp_number_architypes.hpp:

* The file name is misspelled.  It should
  be "archetypes."

* The #includes aren't quite right.  The
  file does not use trunc, and does use
  fpclassify.

* line 28: mpl::list requires #include <boost/mpl/list.hpp>
* line 79: stringstream requires #include <sstream>

depricated:

* Should be deprecated.

detail/functions/constants.hpp:

* line 12: typedef typename has_enough_bits<...> pred_type;
  If you're going to go through all this trouble to make
  sure that you have a type with enough bits, you shouldn't
  assume that unsigned has enough.  IIRC, the standard
  only guarantees 16 bits for unsigned.

* line 38: if(digits  < 3640)
  Ugh.  I'm assuming that digits means binary digits
  and 3640 = floor(log_2 10^1100)

* line 51: I would appreciate a brief description
  of the formula.  This code isn't very readable.

* line 105: It looks like this is using the formula
  e \approx \sum_0^n{\frac{1}{i!}} = (\sum_0^n\frac{n!}{i!})/n!
  and induction using the rule:
  \sum_0^n\frac{n!}{i!} =
    (\sum_0^(n-1)\frac{(n-1)!}{i!}) * n + 1

  Please explain the math or point to somewhere
  that explains it.

* line 147: ditto for pi

* line 251: calc_pi(result, ...digits2<...>::value + 10);
  Any particular reason for +10?  I'd prefer
  calc_pi to do whatever adjustments are
  needed to get the right number of digits.

detail/functions/pow.hpp:

* line 18 ff: pow_imp
  I don't like the fact that this implementation
  creates log2(p) temporaries on the
  stack and multiplies them all together
  on the way out of the recursion.  I'd
  prefer to only use a constant number
  of temporaries.

* line 31:
  I don't see how the cases > 1 buy you anything.
  They don't change the number of calls to eval_multiply.
  (At least, it wouldn't if the algorithm were
  implemented correctly.)

* line 47:
  ???: The work of this loop is repeated at every
  level of recursion. Please fix this function.
  It's a well-known algorithm.  It shouldn't be
  hard to get it right.

* line 76: pow_imp(denom, t, -p, mpl::false_());
  This is broken for p == INT_MIN with twos-complement
  arithmetic.

* line 185: "The ldexp function..."
  I think you mean the "exp" function.

* line 273: if(x.compare(ll) == 0)
  This is only legal if long long is in T::int_types.
  Better to use intmax_t which is guaranteed to be in
  T::int_types.

* line 287: const bool b_scale = (xx.compare(float_type(1e-4)) > 0);
  Unless I'm missing something, this will always be
  true because xx > 1 > 1e-4 as a result of the
  test on line 240.

* line 332: if(&arg == &result)
  This is totally unnecessary, since you already copy
  arg into a temporary using frexp on line 348.

* line 351: if(t.compare(fp_type(2) / fp_type(3)) <= 0)
  Is the rounding error evaluating 2/3 important?
  I've convinced myself that it works, but it
  took a little thought.  Maybe add a comment to
  the effect that the exact value of the boundary
  doesn't matter as long as its about 2/3?

* line 413: "The fabs function "
  s/fabs/log10/.

* line 464: eval_convert_to(&an, a);
  This makes the assumption that conversion
  to long long is always okay even if the
  value is out of the range of long long.  This
  assumption makes me nervous.

* line 465: a.compare(an)
  compare(long long)

* line 474:
  eval_subtract(T, long long)
  Also, this value of da is never used.  It
  gets overwritten on line 480 or 485.

* line 493: (eval_get_sign(x) > 0)
  This expression should always be true at this point.

* line 614: *p_cosh = x;
  cosh is an even function, so this is wrong when
  x = -\inf

* line 625: T e_px, e_mx;
  These are only used in the generic implementation.
  It's better to restrict the scope of variables
  to the region that they're actually used in.

* line 643: eval_divide(*p_sinh, ui_type(2));
  I notice that in other places you use ldexp for
  division by powers of 2.  Which is better?

detail/mp_number_base.hpp:

* line 423: (std::numeric_limits<T>::digits + 1) * 1000L
  This limits the number of bits to LONG_MAX/1000.
  I'd like there to be some way to know that as
  long as I use less than k bits, the library
  won't have any integer overflow problems.  k
  can depend on INT_MAX, LONG_MAX, etc, but the
  implementation limits need to be clearly documented
  somewhere.  I really hate the normal ad hoc
  situation where you just ignore the problem and
  hope that all the values are small enough to avoid
  overflow.

Other notes:

Backend Requirements:

* Do we always require that all number types
  be assignable from intmax_t, uintmax_t,
  and long double?  Why must an integer
  be assignable from long double?

* The code seems to assume that it's okay to
  pass int as the exp parameter of ldexp.
  These requirements don't mandate that this
  will work.

In Christ,
Steven Watanabe