[rnhmjoj at inventati.org: Re: Contributing to GMP]

[Michele Guerini Rocco]

I forgot to add a CC. I'm not used to mailing lists, sorry.

----- Forwarded message from Michele Guerini Rocco <rnhmjoj at inventati.org> -----

Date: Thu, 20 Feb 2020 00:07:06 +0100
From: Michele Guerini Rocco <rnhmjoj at inventati.org>
To: Anders Andersson <pipatron at gmail.com>
Subject: Re: Contributing to GMP

Well, you right, log(N) is mostly surely not going to be a rational:
I used fractions because I was already using them to represent
(exact) rationals in other parts of the program. In fact the precision
of the computation is controlled by the number of decimal places, so it
makes more sense to use floating points.
I guess I could rewrite the function to use floats and the algorithm
should still work; I have not used mpf before, though.
Anyway, here's the description of `mpq_log_ui``


    /* `mpq_log_ui(z, x, D)`
     * calculates the logarithm of x up
     * to D decimal places using the formula
     *
     *   log(x) = log(a(x) * 2^(n-1))
     *          = log(a) + (n-1) log(2)
     *
     * where 1<a<2 and n is the number of binary
     * digits of x.
     *
     * log(a) is the computed from
     *
     *   log((1+y)/(1-y)) = 2Σ_{k=0} y^(2k+1)/(2k+1)
     *
     * where y = (a-1)/(a+1).
     *
     * The error of the n-th partial sum is
     * estimated as:
     *
     *   Δ(n, y) = 4y^(2n + 3) /
     *        (4n^2y^4 - 8n^2y^2 + 4n^2 + 4ny^4 - 12ny^2 + 8n + y^4 - 4y^2 + 3)
     *
     * So to fix the first D decimal places
     * we evaluate up to n terms such that
     *
     *   Δ(n, y) < 1/10^D
     */

log(2) above is computed with


    /* `mpq_log2(z, D)`
     * calculate log(2) to D decimal places
     * using the series
     *
     *   log(2) = Σ_{k=1} 1/(k 2^k)
     *
     * The error of the n-th partial sum is
     * estimated as
     *
     *   Δn = 1/(n(n+2) 2^(n-1)).
     *
     * So to fix the first D decimal places
     * we evaluate up to n terms such that
     *
     *   Δn < 1/10^D
     */

On 18-02-20, Anders Andersson wrote:
> On Tue, Feb 18, 2020 at 1:09 PM Michele Guerini Rocco
> <rnhmjoj at inventati.org> wrote:
> >
> > Hello everyone,
> >
> > I have written, as an exercise, a program to compute the Euler-Mascheroni
> > constant (γ) to arbitrary precision using GMP. It's not meant to beat
> > existing programs (like y-cruncher) or anything like it, but it works.
> > Anyway, I had to compute the natural logarithm of a natural number to
> > arbitrary precision, so I wrote a `mpq_log_ui` function to do that,
> > something which GMP doesn't provide.
> >
> > Since this function could be useful in general, for more pratical
> > purposes that computing a constant, I thought about contributing it
> > back. I'm a layman in numerical analysis and I'm not sure the
> > implementation it's up to the GMP standards but I can work on it.
> > So how can I contribute to GMP?
> 
> Does it really make sense for this to be a rational number function
> instead of an mpf_log_ui? Everything else dealing with rational
> numbers works with *exact* numbers, and it seems to me that a log_ui
> will rarely be exact.


----- End forwarded message -----
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 228 bytes
Desc: not available
URL: <https://gmplib.org/list-archives/gmp-discuss/attachments/20200220/e85700ff/attachment-0002.bin>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 228 bytes
Desc: not available
URL: <https://gmplib.org/list-archives/gmp-discuss/attachments/20200220/e85700ff/attachment-0003.bin>