logo标签列表

有没有GMP对数函数？

gmplogarithm
15

GMP库中是否实现了对数函数？

5个回答
6
没有GMP中的这种函数。 只有在MPFR中。
你是正确的，log函数在gmp库中没有被明确给出，然而通过使用mpz_get_d_2exp并利用对数的性质，可以获得log值。 - Joseph Wood
6
我知道您没有询问如何实现它，但是...
您可以使用对数的属性来粗略地实现它：http://gnumbers.blogspot.com.au/2011/10/logarithm-of-large-number-it-is-not.html 以及 GMP 库的内部：https://gmplib.org/manual/Integer-Internals.html （编辑：基本上只需使用 GMP 表示的最高位“数字”，因为表示的基数是巨大的 B^N 远大于 B^{N-1}）
这是我为 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;
    }

（后来编辑）五年后回来看，我认为很酷的是，在本机浮点实现中甚至出现了log(a) = N log(B) + log(a_N)的核心概念，这里是ia64的glibc实现。在遇到这个问题后，我再次使用了它。
5
正如其他答案所说，GMP中没有对数函数。部分答案提供了对数函数的实现，但仅具有双精度精度，而不是无限精度。
我下面实现了完整（任意）精度的对数函数，甚至可以高达数千位的精度。使用GMP的通用浮点类型mpf。
我的代码使用ln(1 + x)的Taylor系列加上mpf_sqrt()（用于提高计算效率）。
代码使用C++编写，并且由于两个事实而相当庞大。首先，它进行精确的时间测量，以找出您机器的内部计算参数的最佳组合。第二个原因是它使用额外的速度提高，例如利用mpf_sqrt()来准备初始值。
我的代码的算法如下：

  1. 从输入参数 x 中提取 2 的指数幂，即使用 mpf_get_d_2exp() 重写 x = d * 2^exp。

  2. 使得从上一步得到的 d 满足不等式 2/3 <= d <= 4/3，这可以通过将 d 乘以 2 并做 --exp 操作来实现。这确保了 d 和 1 最多相差 1/3，换句话说，d 向正负两个方向都扩展了相同的距离。

  3. 使用 mpf_div_2exp() 和 mpf_mul_2exp() 将 x 除以 2^exp。

  4. 对 x 进行多次（num_sqrt 次）平方根运算，使得 x 接近 1。这将确保 Taylor 级数更快地收敛。因为进行多次平方根运算比在 Taylor 级数的额外迭代中花费更多时间要快。

  5. 计算 ln(1 + x) 的 Taylor 级数 直到所需精度（甚至需要数千位的精度）。

  6. 因为在步骤 4. 中我们多次进行了平方根运算，现在我们需要将 y（Taylor 级数的结果）乘以 2^num_sqrt。

  7. 最后，在步骤 1. 中我们提取了 2^exp，现在我们需要将 ln(2) * exp 加到 y 中。这里的 ln(2) 只需通过同一个实现整个算法的函数进行一次递归调用即可计算出来。

上述步骤源于公式序列：ln(x) = ln(d * 2^exp) = ln(d) + exp * ln(2) = ln(sqrt(...sqrt(d))) * num_sqrt + exp * ln(2)。

我的实现会自动进行定时（每次程序运行仅一次），以确定需要多少次平方根运算来平衡 Taylor 级数的计算。如果您需要避免定时，则可以将 mpf_ln() 的第三个参数 sqrt_range 传递为 0.001，而不是零。

main()函数包含使用示例、正确性测试（通过与低精度 std::log()进行比较）、时间测量和不同详细信息的输出。该函数在π数的前1024位上进行测试。

在调用我的mpf_ln()函数之前，请不要忘记通过调用mpf_set_default_prec(bits)来设置所需的计算精度，以位为单位。

我的mpf_ln()的计算时间约为40-90微秒，适用于1024位精度。更大的精度将需要更多的时间，这大致成正比于精度位数的数量。

函数的第一次运行需要相对较长的时间，因为它执行时序表的预计算和ln(2)值的计算。因此建议在程序启动时首先进行单个计算，以避免后面在代码中出现的时间关键区域内的较长计算。

例如在Linux上编译，您必须安装GMP库并发出以下命令：

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

在线试用!

#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;
}

输出：

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
5
下面的方法利用了 mpz_get_d_2exp，并且是从 gmp R 包 获得的。它可以在文件 bigintegerR.cc 中的函数 biginteger_log 下找到（您首先需要下载源代码（即tar文件））。您还可以在此处查看它：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;
}

当然，上述方法可以进行修改以返回具有任何基数的对数，使用对数的属性（例如，变底公式）。
2
优秀的解决方案！ - Arty
3

这是它的链接： https://github.com/linas/anant

提供 GNU MP 实数和复数对数、指数、正弦、余弦、伽马函数、反正切函数、平方根、多对数 Riemann 和 Hurwitz zeta 函数、互变超几何函数、拓扑正弦函数等更多功能。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接