Arb 任意精度计算库

             

预备知识 C++ 基础

   Arb 是一款支持任意精度计算的 C/C++ 程序,支持对计算进行严谨的误差估计,即每个任意精度浮点数 $z$ 都会伴随一个误差半径 $r$,使得精确结果必定落在复平面上以 $z$ 为圆心半径为 $r$ 的圆盘中.通过增加浮点数的精度,就可以用数值方法无限逼近理论值.Arb 还提供了许多特殊函数的计算,例如 $\Gamma$ 函数,超几何函数 等,以及常用的线性代数功能和离散傅里叶变换等.官方主页 arblib.org 包含详细的文档.

1. 安装

   以 Ubuntu 为例,最方便的安装方式就是使用 apt 安装.首先安装 dependency

sudo apt install libmpfr-dev libflint-dev
其中 MPFRflint 两个包分别用于任意精度浮点数以及数论.然后安装 Arb
sudo apt install libflint-arb-dev libflint-arb2

   但目前这并不是 Arb 的最新版本(例如没有实现库仑函数).要获得最新版本,可以直接从 GitHub 下载源码编译即可(默认使用 gcc 编译器).目前笔者使用的版本是 release 2.19.0.

   可以用 ./configure --help 查看编译选项,若所有的包都安装在默认目录则不需要编译选项.

./configure [编译选项];
make -j4;
sudo make install;
其中 -j4 是使用 4 线程进行编译,也可以改成其他数字.

2. 编译

   在 Ubuntu 中如果你使用 apt 安装,在编译程序是需要加 -lflint-arb 选项.如果你直接从源码编译,则需要加 -larb 选项.

   一个 C++ 例程:该程序用 80 bit 初始精度计算复参数的超几何函数 $_1F_1$ 并于 Mathematica 的结果比较.如果 Arb 估计的结果精度小于 15 位有效数字则显示警告.初始精度越高,结果的有效数字也越高,具体取决于 $_1F_1$ 的参数.读者可以尝试在程序内加入一个循环,若结果有效数字不够,则把初始精度加倍再次计算,直到达到精度为止.

代码 1:test.cpp
#include "arb_hypgeom.h"
#include "acb_hypgeom.h"
#include <complex>
#include <iostream>

using namespace std;

typedef double Doub, Doub_O;
typedef const double Doub_I;
typedef complex<double> Comp, Comp_O;
typedef complex<double> Comp_I;
typedef int Int;
typedef const int Int_I;
typedef int & Int_O, Int_IO;

// 1F1 hypdergeometric function with complex arguments
Comp hypergeom1F1(Comp_I a, Comp_I b, Comp_I z)
{
	slong prec = 80; // set precision bit (log10/log2 = 3.322)
	Comp res;
	arb_t temp1;
	arb_init(temp1);
	acb_t a1, b1, z1, res1;
	acb_init(a1); acb_init(b1); acb_init(z1); acb_init(res1);
	acb_set_d_d(a1, real(a), imag(a));
	acb_set_d_d(b1, real(b), imag(b));
	acb_set_d_d(z1, real(z), imag(z));
	acb_hypgeom_1f1(res1, a1, b1, z1, 0, prec);
	// acb_printn(res1, 100, 0); printf("\n"); // print result
	int digits = acb_rel_accuracy_bits(res1)/3.321928;
	if (digits < 15) {
		cout << "warning: hypergeom1F1 error too large" << endl;
	}
	acb_get_real(temp1, res1);
	res.real(arf_get_d(arb_midref(temp1), ARF_RND_NEAR));
	acb_get_imag(temp1, res1);
	res.imag(arf_get_d(arb_midref(temp1), ARF_RND_NEAR));
	acb_clear(a1); acb_clear(b1); acb_clear(z1);
    acb_clear(res1); arb_clear(temp1);
	return res;
}

int main()
{
	cout << "hypergeom1F1(1.23+1.23I, 1.23+1.23I, 1.23+1.23I) = " << endl;
	cout << hypergeom1F1(Comp(1.23,1.23),Comp(1.23,1.23),Comp(1.23,1.23))
    cout << endl;
	printf("Mathematica: 1.143503984180676 + 3.224470526790991i\n");
}
编译:
g++ -o test.x test.cpp -larb -lflint

致读者: 小时百科一直以来坚持所有内容免费无广告,这导致我们处于严重的亏损状态。 长此以往很可能会最终导致我们不得不选择会员制,大量广告,内容付费等。 因此,我们请求广大读者热心打赏 ,使网站得以健康发展。 如果看到这条信息的每位读者能慷慨打赏 10 元,我们一个星期内就能脱离亏损, 并保证网站能在接下来的一整年里向所有读者继续免费提供优质内容。 但遗憾的是只有不到 1% 的读者愿意捐款, 他们的付出帮助了 99% 的读者免费获取知识, 我们在此表示感谢。

广告位

投放详情

         

© 小时科技 保留一切权利