Is there any GMP logarithm function?

Question

Is there any logarithm function implemented in the GMP library?

Answer 1

The method below makes use of mpz_get_d_2exp and was obtained from the gmp R package. It can be found under the function biginteger_log in the file bigintegerR.cc (You first have to download the source (i.e. the tar file)). You can also see it here: biginteger_log.

// Adapted for general use from the original biginteger_log
// xi = di * 2 ^ ex  ==> log(xi) = log(di) + ex * log(2)

double biginteger_log_modified(mpz_t x) {
    signed long int ex;
    const double di = mpz_get_d_2exp(&ex, x);
    return log(di) + log(2) * (double) ex;
}

Of course, the above method could be modified to return the log with any base using the properties of logarithm (e.g. the change of base formula).

Answer 2

As other answers said, there is no logarithmic function in GMP. Part of answers provided implementations of logarithmic function, but with double precision only, not infinite precision.

I implemented full (arbitrary) precision logarithmic function below, even up to thousands bits of precision if you wish. Using mpf, generic floating point type of GMP.

My code uses Taylor serie for ln(1 + x) plus mpf_sqrt() (for boosting computation).

Code is in C++, and is quite large due to two facts. First is that it does precise time measurements to figure out best combinations of internal computational parameters for your machine. Second is that it uses extra speed improvements like extra usage of mpf_sqrt() for preparing initial value.

Algorithm of my code is following:

1. Factor out exponent of 2 from input x, i.e. rewrite x = d * 2^exp, with usage of mpf_get_d_2exp().

2. Make d (from step above) such that 2/3 <= d <= 4/3, this is achieved by possibly multiplying d by 2 and doing --exp. This ensures that d always differs from 1 by at most 1/3, in other words d extends from 1 in both directions (negative and positive) in equal distance.

3. Divide x by 2^exp, with usage of mpf_div_2exp() and mpf_mul_2exp().

4. Take square root of x several times (num_sqrt times) so that x becomes closer to 1. This ensures that Taylor Serie converges more rapidly. Because computation of square root several times is faster than contributing much more time in extra iterations of Taylor Serie.

5. Compute Taylor Serie for ln(1 + x) up to desired precision (even thousands of bit of precision if needed).

6. Because in Step 4. we took square root several times, now we need to multiply y (result of Taylor Serie) by 2^num_sqrt.

7. Finally because in Step 1. we factored out 2^exp, now we need to add ln(2) * exp to y. Here ln(2) is computed by just one recursive call to same function that implements whole algorithm.

Steps above come from sequence of formulas ln(x) = ln(d * 2^exp) = ln(d) + exp * ln(2) = ln(sqrt(...sqrt(d))) * num_sqrt + exp * ln(2).

My implementation automatically does timings (just once per program run) to figure out how many square roots is needed to balance out Taylor Serie computation. If you need to avoid timings then pass 3rd parameter sqrt_range to mpf_ln() equal to 0.001 instead of zero.

main() function contains examples of usage, testing of correctness (by comparing to lower precision std::log()), timings and output of different verbose information. Function is tested on first 1024 bits of Pi number.

Before call to my function mpf_ln() don't forget to setup needed precision of computation by calling mpf_set_default_prec(bits) with desired precision in bits.

Computational time of my mpf_ln() is about 40-90 micro-seconds for 1024 bit precision. Bigger precision will take more time, that is approximately linearly proportional to the amount of precision bits.

Very first run of a function takes considerably longer time becuse it does pre-computation of timings table and value of ln(2). So it is suggested to do first single computation at program start to avoid longer computation inside time critical region later in code.

To compile for example on Linux, you have to install GMP library and issue command:

clang++-14 -std=c++20 -O3 -lgmp -lgmpxx -o main main.cpp && ./main

Try it online!

#include <cstdint>
#include <iomanip>
#include <iostream>
#include <cmath>
#include <chrono>
#include <mutex>
#include <vector>
#include <unordered_map>

#include <gmpxx.h>

double Time() {
    static auto const gtb = std::chrono::high_resolution_clock::now();
    return std::chrono::duration_cast<std::chrono::duration<double>>(
        std::chrono::high_resolution_clock::now() - gtb).count();
}

mpf_class mpf_ln(mpf_class x, bool verbose = false, double sqrt_range = 0) {
    auto total_time = verbose ? Time() : 0.0;
    int const prec = mpf_get_prec(x.get_mpf_t());
    
    if (sqrt_range == 0) {
        static std::mutex mux;
        std::lock_guard<std::mutex> lock(mux);
        static std::vector<std::pair<size_t, double>> ranges;
        if (ranges.empty())
            mpf_ln(3.14, false, 0.01);
        while (ranges.empty() || ranges.back().first < prec) {
            size_t const bits = ranges.empty() ? 64 : ranges.back().first * 3 / 2;
            mpf_class x = 3.14;
            mpf_set_prec(x.get_mpf_t(), bits);
            double sr = 0.35, sr_best = 1, time_best = 1000;
            size_t constexpr ntests = 5;
            while (true) {
                auto tim = Time();
                for (size_t i = 0; i < ntests; ++i)
                    mpf_ln(x, false, sr);
                tim = (Time() - tim) / ntests;
                bool updated = false;
                if (tim < time_best) {
                    sr_best = sr;
                    time_best = tim;
                    updated = true;
                }
                sr /= 1.5;
                if (sr <= 1e-8) {
                    ranges.push_back(std::make_pair(bits, sr_best));
                    break;
                }
            }
        }
        sqrt_range = std::lower_bound(ranges.begin(), ranges.end(), size_t(prec),
            [](auto const & a, auto const & b){
                return a.first < b;
            })->second;
    }

    signed long int exp = 0;
    // https://gmplib.org/manual/Converting-Floats
    double d = mpf_get_d_2exp(&exp, x.get_mpf_t());
    if (d < 2.0 / 3) {
        d *= 2;
        --exp;
    }
    mpf_class t;
    // https://gmplib.org/manual/Float-Arithmetic
    if (exp >= 0)
        mpf_div_2exp(x.get_mpf_t(), x.get_mpf_t(), exp);
    else
        mpf_mul_2exp(x.get_mpf_t(), x.get_mpf_t(), -exp);
    auto sqrt_time = verbose ? Time() : 0.0;
    // Multiple Sqrt of x
    int num_sqrt = 0;
    if (x >= 1)
        while (x >= 1.0 + sqrt_range) {
            // https://gmplib.org/manual/Float-Arithmetic
            mpf_sqrt(x.get_mpf_t(), x.get_mpf_t());
            ++num_sqrt;
        }
    else
        while (x <= 1.0 - sqrt_range) {
            mpf_sqrt(x.get_mpf_t(), x.get_mpf_t());
            ++num_sqrt;
        }
    if (verbose)
        sqrt_time = Time() - sqrt_time;

    static mpf_class const eps = [&]{
        mpf_class eps = 1;
        mpf_div_2exp(eps.get_mpf_t(), eps.get_mpf_t(), prec + 8);
        return eps;
    }(), meps = -eps;

    // Taylor Serie for ln(1 + x)
    // https://math.stackexchange.com/a/878376/826258
    x -= 1;
    mpf_class k = x, y = x, mx = -x;
    size_t num_iters = 0;
    for (int32_t i = 2;; ++i) {
        k *= mx;
        y += k / i;
        // Check if error is small enough
        if (meps <= k && k <= eps) {
            num_iters = i;
            break;
        }
    }
    auto VerboseInfo = [&]{
        if (!verbose)
            return;
        total_time = Time() - total_time;
        std::cout << std::fixed << "Sqrt range " << sqrt_range << ", num sqrts "
            << num_sqrt << ", sqrt time " << sqrt_time << " sec" << std::endl;
        std::cout << "Ln number of iterations " << num_iters << ", ln time "
            << total_time << " sec" << std::endl;
    };
    // Correction due to multiple sqrt of x
    y *= 1 << num_sqrt;
    if (exp == 0) {
        VerboseInfo();
        return y;
    }
    mpf_class ln2;
    {
        static std::mutex mutex;
        std::lock_guard<std::mutex> lock(mutex);
        static std::unordered_map<size_t, mpf_class> ln2s;
        auto it = ln2s.find(size_t(prec));
        if (it == ln2s.end()) {
            mpf_class sqrt_sqrt_2 = 2;
            mpf_sqrt(sqrt_sqrt_2.get_mpf_t(), sqrt_sqrt_2.get_mpf_t());
            mpf_sqrt(sqrt_sqrt_2.get_mpf_t(), sqrt_sqrt_2.get_mpf_t());
            it = ln2s.insert(std::make_pair(size_t(prec), mpf_class(mpf_ln(sqrt_sqrt_2, false, sqrt_range) * 4))).first;
        }
        ln2 = it->second;
    }
    y += ln2 * exp;
    VerboseInfo();
    return y;
}

std::string mpf_str(mpf_class const & x) {
    mp_exp_t exp;
    auto s = x.get_str(exp);
    return s.substr(0, exp) + "." + s.substr(exp);
}

int main() {
    // https://gmplib.org/manual/Initializing-Floats
    mpf_set_default_prec(1024); // bit-precision
    // http://www.math.com/tables/constants/pi.htm
    mpf_class x(
        "3."
        "1415926535 8979323846 2643383279 5028841971 6939937510 "
        "5820974944 5923078164 0628620899 8628034825 3421170679 "
        "8214808651 3282306647 0938446095 5058223172 5359408128 "
        "4811174502 8410270193 8521105559 6446229489 5493038196 "
        "4428810975 6659334461 2847564823 3786783165 2712019091 "
        "4564856692 3460348610 4543266482 1339360726 0249141273 "
        "7245870066 0631558817 4881520920 9628292540 9171536436 "
    );
    std::cout << std::boolalpha << std::fixed << std::setprecision(14);
    std::cout << "x:" << std::endl << mpf_str(x) << std::endl;
    auto cmath_val = std::log(mpf_get_d(x.get_mpf_t()));
    std::cout << "cmath ln(x): " << std::endl << cmath_val << std::endl;
    auto volatile tmp = mpf_ln(x); // Pre-Compute to heat-up timings table.
    auto time_start = Time();
    size_t constexpr ntests = 20;
    for (size_t i = 0; i < ntests; ++i) {
        auto volatile tmp = mpf_ln(x);
    }
    std::cout << "mpf ln(x) time " << (Time() - time_start) / ntests << " sec" << std::endl;
    auto mpf_val = mpf_ln(x, true);
    std::cout << "mpf ln(x):" << std::endl << mpf_str(mpf_val) << std::endl;
    std::cout << "equal to cmath: " << (std::abs(mpf_get_d(mpf_val.get_mpf_t()) - cmath_val) <= 1e-14) << std::endl;
    return 0;
}

Output:

x:
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587007
cmath ln(x): 
1.14472988584940
mpf ln(x) time 0.00004426845000 sec
Sqrt range 0.00000004747981, num sqrts 23, sqrt time 0.00001440000000 sec
Ln number of iterations 42, ln time 0.00003873100000 sec
mpf ln(x):
1.144729885849400174143427351353058711647294812915311571513623071472137769884826079783623270275489707702009812228697989159048205527923456587279081078810286825276393914266345902902484773358869937789203119630824756794011916028217227379888126563178049823697313310695003600064405487263880223270096433504959511813198
equal to cmath: true

Answer 3

I know you didn't ask how to implement it, but...

You can implement a rough one using the properties of logarithms: http://gnumbers.blogspot.com.au/2011/10/logarithm-of-large-number-it-is-not.html

And the internals of the GMP library: https://gmplib.org/manual/Integer-Internals.html

(Edit: Basically you just use the most significant "digit" of the GMP representation since the base of the representation is huge B^N is much larger than B^{N-1})

Here is my implementation for Rationals.

    double LogE(mpq_t m_op)
    {
        // log(a/b) = log(a) - log(b)
        // And if a is represented in base B as:
        // a = a_N B^N + a_{N-1} B^{N-1} + ... + a_0
        // => log(a) \approx log(a_N B^N)
        // = log(a_N) + N log(B)
        // where B is the base; ie: ULONG_MAX

        static double logB = log(ULONG_MAX);

        // Undefined logs (should probably return NAN in second case?)
        if (mpz_get_ui(mpq_numref(m_op)) == 0 || mpz_sgn(mpq_numref(m_op)) < 0)
            return -INFINITY;               

        // Log of numerator
        double lognum = log(mpq_numref(m_op)->_mp_d[abs(mpq_numref(m_op)->_mp_size) - 1]);
        lognum += (abs(mpq_numref(m_op)->_mp_size)-1) * logB;

        // Subtract log of denominator, if it exists
        if (abs(mpq_denref(m_op)->_mp_size) > 0)
        {
            lognum -= log(mpq_denref(m_op)->_mp_d[abs(mpq_denref(m_op)->_mp_size)-1]);
            lognum -= (abs(mpq_denref(m_op)->_mp_size)-1) * logB;
        }
        return lognum;
    }

(Much later edit) Coming back to this 5 years later, I just think it's cool that the core concept of log(a) = N log(B) + log(a_N) shows up even in native floating point implementations, here is the glibc one for ia64 And I used it again after encountering this question

Answer 4

No there is no such function in GMP. Only in MPFR.

Answer 5

Here it is: https://github.com/linas/anant

Provides gnu mp real and complex logarithm, exp, sine, cosine, gamma, arctan, sqrt, polylogarithm Riemann and Hurwitz zeta, confluent hypergeometric, topologists sine, and more.