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

Subject: Re: [boost] [review] Multiprecision review (June 8th - 17th, 2012)
From: John Maddock (boost.regex_at_[hidden])
Date: 2012-06-25 08:19:53

> 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.

Many thanks Steven.

> concepts/mp_number_architypes.hpp:
>
> * The file name is misspelled.  It should
>  be "archetypes."

Indeed, fixed in my local copy.

> * 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>

Fixed in my local copy.

> depricated:
>
> * Should be deprecated.

Will go away if the library is accepted anyway (just contains bits and bobs 
not in the submission that I didn't want to delete just yet).

> 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.

Nod.  There's some new code that simplifies all that as well, simplified and 
fixed in my local copy.

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

Comments to that effect throughout.

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

Nod.  Added locally:

   //
   // We calculate log2 from using the formula:
   //
   // ln(2) = 3/4 SUM[n>=0] ((-1)^n * N!^2 / (2^n(2n+1)!))
   //
   // Numerator and denominator are calculated separately and then
   // divided at the end, we also precalculate the terms up to n = 5
   // since these fit in a 32-bit integer anyway.
   //
   // See Gourdon, X., and Sebah, P. The logarithmic constant: log 2, Jan. 
2004.
   // Also http://www.mpfr.org/algorithms.pdf.
   //

Also spotted a missing optimisation in the arithmetic.

> * 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.

Nod adding locally:

   //
   // Standard evaluation from the definition of e: 
http://functions.wolfram.com/Constants/E/02/
   //

> * line 147: ditto for pi

Nod.  Adding locally:

   //
   // This algorithm is from:
   // Schonhage, A., Grotefeld, A. F. W., and Vetter, E. Fast Algorithms: A 
Multitape Turing
   // Machine Implementation. BI Wissenschaftverlag, 1994.
   // Also described in MPFR's algorithm guide: 
http://www.mpfr.org/algorithms.pdf.
   //
   // Let:
   // a[0] = A[0] = 1
   // B[0] = 1/2
   // D[0] = 1/4
   // Then:
   // S[k+1] = (A[k]+B[k]) / 4
   // b[k] = sqrt(B[k])
   // a[k+1] = a[k]^2
   // B[k+1] = 2(A[k+1]-S[k+1])
   // D[k+1] = D[k] - 2^k(A[k+1]-B[k+1])
   // Stop when |A[k]-B[k]| <= 2^(k-p)
   // and PI = B[k]/D[k]

Currently testing those changes before moving on to the other comments,

Thanks again, John.