开发者

gmp pow with two mpf_t

Is there an implementation in gmp that allows a power function with only mpf_t's as argument? I want to do this:

mpf_t s ;
mpf_init (s);
mpf_set_d (s,boost::lexical_cast<double>(sec));
mpf_t ten,mil;
mpf_init(ten);
mpf_ini开发者_如何学Got(mil);
mpf_set_d(ten,10.0);
mpf_set_d(mil,0.001);
mpf_div(s,s,ten);
mpf_pow_ui(s,ten,s); //<- this doesn't work because it need an unsigned int as third argument but I need it with a mpf_t
mpf_mul(s,s,mil);


I don't think so, at least not with GNU Multi-Precision library only. But you could use mpfr, which is based on gmp and supports a mpfr_pow (mpfr_t rop, mpfr_t op1, mpfr_t op2, mpfr_rnd_t rnd) function. See here.

If you decide to do that, this could also be helpful.


There is one interesting workaround using square root mpf_sqrt_ui. From math we know that x^y = Sqrt(x)^(y * 2), so we can multiply Y many times by 2 and take square root of X same amount of times.

Thus by multiplying Y by 2 you may make it almost whole integer. And as you know there is mpf_pow_ui that does powering into whole integer.

Following code does all this. Don't forget that b should be set to high precision only to allow many times square rooting.

For simplicity I used mpf_class, this is C++ interface to mpf.

I did output to console actual mpf result value and reference value computed through std::pow from .

To avoid setting high precision is a bit more difficult, but possible e.g. through Taylor Serie like

Sqrt(1 + x) = 1 + 1/2*x - 1/8*x^2 + 1/16*x^3 - 5/128*x^4 + ...

Try it online!

#include <cmath>
#include <iostream>
#include <iomanip>

#include <gmpxx.h>

int main() {
    mpf_class const b0 = 9.87654321, p0 = 1.23456789;
    mpf_class b = b0, p = p0;
    b.set_prec(1 << 7); p.set_prec(1 << 7);
    int const sqrt_cnt = 48;
    for (int i = 0; i < sqrt_cnt; ++i)
        mpf_sqrt(b.get_mpf_t(), b.get_mpf_t());
    mpf_mul_2exp(p.get_mpf_t(), p.get_mpf_t(), sqrt_cnt);
    mpf_pow_ui(b.get_mpf_t(), b.get_mpf_t(), std::lround(p.get_d()));
    std::cout << std::fixed << std::setprecision(12) << "Actual "
        << b.get_d() << ", Reference " << std::pow(b0.get_d(), p0.get_d())
        << std::endl;
}

Output:

Actual 16.900803674719, Reference 16.900803674719
0

上一篇:

下一篇:

精彩评论

暂无评论...
验证码 换一张
取 消

最新问答

问答排行榜